TP MP* n°5 (SCILAB) Discrétisation d'équations différentielles

funcprot(0)
// euler_exp impératif
function  E = euler_imp(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for k = 1:N
    un = un + h*f(tn,un)
    tn = tn + h
    X = [X,tn]
    Y = [Y,un]
  end
  E = [X;Y]
endfunction
 
 
function plot_euler_imp(f,a,b,N,yo)
  e = euler_imp(f,a,b,N,yo)
  plot(e(1,:) , e(2,:))
endfunction
 
// euler_exp récursif
function  E = euler_rec_les_y(f,x,y,xmax,h,liste)
  if x > xmax then E = liste
  else E = euler_rec_les_y(f,x+h,y+h*f(x,y),xmax,h,[e,y+h*f(x,y)])
  end
endfunction
 
function euler_rec(f,a,b,N,yo)
  h = (b-a)/N
  x = [a:h:b]
  plot(x,euler_rec_les_y(f,a,yo,b,h,[yo]))
endfunction
 
//plot_euler_imp(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")
 
//euler_rec(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")
 
 
//max(abs(euler_rec_les_y(deff('[u] = f(x,y)', 'u = y'),0,1,2,0.01,[1]) - exp([0:0.01:2])))
 
//e = euler_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//max(abs((e(2,:) - exp(e(1,:)))))
 
 
 
// euler2_exp impératif
function  E = euler2_imp(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for k = 1:N do
    un = un + h*f(tn + 0.5*h , un + 0.5*h*f(tn,un) )
    tn = tn + h
    X = [X,tn]
    Y = [Y,un]
  end
  E = [X;Y]
endfunction
 
function plot_euler2_imp(f,a,b,N,yo)
  e = euler2_imp(f,a,b,N,yo)
  plot(e(1,:) , e(2,:))
endfunction
 
// euler2_exp récursif
function  E = euler2_rec_les_y(f,x,y,xmax,h,liste)
  if x > xmax then E = liste
  else E = euler2_rec_les_y(f,x+h,y+h*f(x+0.5*h,y+0.5*h*f(x,y)),xmax,h,[e,y+h*f(x+0.5*h,y+0.5*h*f(x,y))])
  end
endfunction
 
function euler2_rec(f,a,b,N,yo)
  h = (b-a)/N
  x = [a:h:b]
  plot(x,euler2_rec_les_y(f,a,yo,b,h,[yo]))
endfunction
 
 
//plot_euler2_imp(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")
 
//euler2_rec(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")
 
 
//max(abs(euler2_rec_les_y(deff('[u] = f(x,y)', 'u = y'),0,1,2,0.01,[1]) - exp([0:0.01:2])))
//e1 = euler_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//e2 = euler2_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//max(abs((e1(2,:) - exp(e1(1,:)))))
//max(abs((e2(2,:) - exp(e2(1,:)))))
 
 
//RK4
 
 
function R = rk4(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for i = 1:N
    k1 = h*f(tn , un)
    k2 = h*f(tn+0.5*h , un+k1/2)
    k3 = h*f(tn+0.5*h , un+k2/2)
    k4 = h*f(tn+h , un+k3)
    tn = tn + h
    un = un + (k1 + 2*k2 + 2*k3 + k4)/6
    X = [X,tn]
    Y = [Y,un]
  end
  R = [X;Y]
endfunction
 
//deff('[f] = f(x,y)', 'f = -150*y + 30');
 
//e = euler_imp(f,0,1,50,1/5+1D-15);
//e2 = euler2_imp(f,0,1,50,1/5+1D-15);
//r = rk4(f,0,1,50,1/5+1D-15);
//plot(r(1,:),r(2,:))
//max(abs((r(2,:) - exp(r(1,:)))))
 
 
//t = 0:0.02:1;
//uo = 1/5+1D-15;
//to = 0;
//u = ode(uo,to,t,f);
//plot(t,u)
 
 
// systèmes et ordre 2
// y1' = 3y1 -4y2
// y'2 = y1 + 3y2
// y1(0) = 2  y2(0) = 0
 
 
//deff('dy=f(t,y)' , 'dy(1) = 3*y(1) - 4*y(2) ; dy(2) = y(1) + 3*y(2)');
//t0 = 0 ; y0 = [2 ; 0] ; t = 0: 0.1 : 2 ;
//z = ode(y0,t0,t,f) ;
//plot(t,z(1,:),t,z(2,:))
//plot(t,2*cos(2*t) .* exp(3*t),"ro")
 
//pendule
//deff('dy=f(t,y)' , 'dy(1) = y(2) ; dy(2) = -y(2) -sin(y(1))');
//t0 = 0 ; y0 = [0 ; 1] ; t = 0: 0.1 : 10 ;
 
//z = ode(y0,t0,t,f) ;
//plot(t,z(1,:),t,z(2,:))
 
 
// portrait de phase
//plot(z(1,:),z(2,:))
 
 
 
// stiff
deff('[f] = f(x,y)', 'f = -1000*y +1000*exp(-x) - exp(-x)');
 
e = euler_imp(f,0,1.5,10,0);
e2 = euler2_imp(f,0,1.5,20,0);
r = rk4(f,0,1.5,20,0);
plot(e(1,:),e(2,:))
 
t = 0:0.005:1.5;
uo = 0;
to = 0;
u = ode(uo,to,t,f);
//plot(t,u)