Méthode d'Euler - Equations du 1er ordre
%-----------------------------------------------------
%% Exemple 1 : $(2y) y' = 1-t $
%-----------------------------------------------------
clear
clc
z = dsolve('2*y*Dy=1-t','y(0)=1');
z = simple(z(1))
n=10
a=0;
b=3;
y0=1
dx=(b-a)/n;
x(1) = a;
ye(1) = y0; % explicite
yi(1) = y0; % implicite
ys(1) = y0; % second ordre
fe = @(x,y) y + dx*(1-x)/(2*y);
fi = @(x,y) 0.5*(y + sqrt(y^2 + 2*dx*(1-x)));
fs = @(x,y) fe(x,y)-0.5*(dx^2)*( 1/(2*y) + ((1-x)^2)/(4*y^3));
for i=1:n
x(i+1) = x(i) + dx;
ye(i+1) = fe(x(i),ye(i));
yi(i+1) = fi(x(i+1),yi(i));
ys(i+1) = fs(x(i),ys(i));
end
ezplot(z(1),[a,b])
hold on
plot(x,ye,'k -.')
plot(x,yi,'k --')
plot(x,ys,'k -.o')
plot(x,0.5*(yi+ye),'k -.+')
Méthode d'Euler - Système de 2 équations
clearclc
a=0;
b=30;
y0=1
dy0=-1
f = inline('[y(2);-.2*y(2)-sin(y(1))+sin(t)]','t','y');
[T Z] = ode45(f,[a b],[y0;dy0]);
%plot(T,Z(:,1),'r',T,Z(:,2))
plot(T,Z(:,1),'r')
%%
n=500;
dx=(b-a)/n;
x(1) = a;
Y = [y0 dy0];
Fe = @(t,V) V + dx*[ V(2);
-.2*V(2)-sin(V(1))+sin(t)];
for i=1:n
x(i+1) = x(i) + dx;
Y(i+1,:) = Fe(x(i),Y(i,:)');
end
plot(x,Y(:,1),'k -.',T,Z(:,1))
Méthode d'accélération moyenne constante (MAMC)
%---------------------------------------
function Ut = MAMC(M,C,K,F,Acc,dt)
%---------------------------------------
K = K + 4*M/dt^2 + 2*C/dt ;
C = 2*C + 4*M/dt;
N = length(Acc);
U = full(0*F);
V = U;
A = U;
Ut(:,1) = U;
for i = 1:N-1
dP = F * (Acc(i+1)-Acc(i));
dP = dP + C * V + 2*M * A;
dU = K\dP;
dA = 4*dU/dt^2 - 4*V/dt - 2*A;
C = 2*C + 4*M/dt;
N = length(Acc);
U = full(0*F);
V = U;
A = U;
Ut(:,1) = U;
for i = 1:N-1
dP = F * (Acc(i+1)-Acc(i));
dP = dP + C * V + 2*M * A;
dU = K\dP;
dA = 4*dU/dt^2 - 4*V/dt - 2*A;
dV = 2*dU/dt - 2*V;
U = U + dU;
V = V + dV;
A = A + dA;
Ut(:,i+1) = U;
end
return
