%============ Resolution numerique des equations differentielles =========
clear all;close all;clc;
a=input ('a =');
b=input ('b =');
n=input ('n =');
h = (b-a)/n;
t=a:h:b;
disp(['h=',num2str(h)])
disp(['*Resolution numerique de l''equation differentielle f''=y-(2.x/y)*'])
disp(['*** S=1. Pour la methode d''Euler ***'])
disp(['*** S=2. Pour la methode de Runge-kutta ***'])
disp(['*** S=3. Pour les deux methodes ***'])
disp(['*** S=4. Pour Quitter ***'])
disp(['*** Tapez S entre ˂1-4˃ ***'])
s=input('S =');
switch s
case 1
disp(['*** Debut de la methode de Runge-Kutta ***'])
%%%%%%%%%%%%%%%%%%%%%%%%%% Trace de la solution exacte %%%%%%%%%%%%%
yex=sqrt(2*t+1);
figure('color',[1 1 1])
plot(t,yex,'*-g')
%%%%%%%%%%%%%%%%%%%%%%%%%%% METHODE Runge-Kutta %%%%%%%%%%%%%%%%%%%%%%%%%
dydt=inline('y-(2.*t/y)','t','y');
epsilon=0.0001;u(1)=1+epsilon;
for i=1:n-1
u1(i)=u(i);u2(i)=u(i)+h/2*dydt(t(i),u1(i));
u3(i)=u(i)+h/2*dydt(t(i)+h/2,u2(i));u4(i)=u(i)+h*dydt(t(i)+h/2,u3(i));
u(i+1)=u(i)+h/6*(dydt(t(i),u1(i))+2*dydt(t(i)+h/2,...
u2(i))+2*dydt(t(i)+h/2,u3(i))+dydt(t(i+1),u4(i)));
end
hold on
plot(t(1:end-1),u,'*-r')
xlabel(' tn')
ylabel(' Un')
legend('Exacte', 'Runge-Kutta')
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%% METHODE d'Euler %%%%%%%%%%%%%%%%%%%%%%%%%
case 2
disp(['*** Debut de la methode d''Euler ***'])
%%%%%%%%%%%%%%%%%%%%%%%%%% Trace de la solution exacte %%%%%%%%%%%%%
yex=sqrt(2*t+1);
figure('color',[1 1 1])
plot(t,yex,'*-g')
dydt=inline('y-(2.*t/y)','t','y');
epsilon=0.0001;u(1)=1+epsilon;
for i=1:n-1
u(i+1)=u(i)+h/2*(dydt(t(i),u(i)));
end
hold on
plot(t(1:end-1),u,'*-b')
xlabel(' tn')
ylabel(' Un')
legend('Exacte', 'Euler')
grid
case 3
disp(['*** Pour les deux methodes ***'])
%%%%%%%%%%%%%%%%%%%%%%%%%% Trace de la solution exacte %%%%%%%%%%%%%
yex=sqrt(2*t+1);
figure('color',[1 1 1])
plot(t,yex,'*-g')
%%%%%%%%%%%%%%%%%%%%%%%%%%% METHODE Runge-Kutta %%%%%%%%%%%%%%%%%%%%%%%%%
dydt=inline('y-(2.*t/y)','t','y');
epsilon=0.0001;u(1)=1+epsilon;
for i=1:n-1
u1(i)=u(i);u2(i)=u(i)+h/2*dydt(t(i),u1(i));
u3(i)=u(i)+h/2*dydt(t(i)+h/2,u2(i));u4(i)=u(i)+h*dydt(t(i)+h/2,u3(i));
u(i+1)=u(i)+h/6*(dydt(t(i),u1(i))+2*dydt(t(i)+h/2,...
u2(i))+2*dydt(t(i)+h/2,u3(i))+dydt(t(i+1),u4(i)));
end
hold on
plot(t(1:end-1),u,'*-r')
%%%%%%%%%%%%%%%%%%%%%%%%%%% METHODE d'Euler %%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n-1
u(i+1)=u(i)+h/2*(dydt(t(i),u(i)));
end
hold on
plot(t(1:end-1),u,'*-b')
xlabel(' tn')
ylabel(' Un')
legend('Exacte', 'Runge-Kutta', 'Euler')
grid
case 4
quit
otherwise
disp('Erreur! la valeur de S entre ˂1-4˃');
end