assignment 3 2.pdf
Bonjour a tous,
Merci de regarder le pdf ci-joint, et on peut echanger des idées dans ce sujet. Merci d'avance.
> restart;
> f := (x, y) ->1/(1+sin(y)) ;
> ode := diff(y(x), x) = f(x, y); x \in [0,x0]
## Algorithme de Huen pour un pas h= x0/(2N); on met 2N points sur le segment [0,X0].
Notre but est de trouver l'erreur ( voir svp) sa relation dans le fichier pdf attaché.
## Methode de Heun avec un pas 2h cad y_i^(2h) aux points x_i=i(2h)
> Heun1 := proc (f, x0, h)
local x, y, i, N, k;
N := round((1/2)*x0/h);
y := Array(0 .. N);
x := Array(0 .. N);
x[0] := 0;
y[0] := (1/4)*Pi;
for i from 0 to N-1 do
x[i+1] := (2*i+2)*h;
k[1] := f(x[i], y[i]);
k[2] := f(x[i]+2*h, y[i]+2*h*k[1]);
y[i+1] := y[i]+h*(k[1]+k[2])
end do;
evalf([seq([x[i], y[i]], i = 0 .. N)])
end proc;
## Methode de Heun avec un pas h cad y_2i^(h) mais aux noeuds x_{2i}=2i*h
Heun2 := proc (f, x0, h)
local x, y, i, N, k;
N := round((1/2)*x0/h);
y := Array(0 .. 2*N);
x := Array(0 .. 2*N);
x[0] := 0;
y[0] := (1/4)*Pi;
for i from 0 to N-1 do
x[2*i+2] := (2*i+2)*h;
k[1] := f(x[2*i], y[2*i]);
k[2] := f(x[2*i]+h, y[2*i]+h*k[1]);
y[2*i+2] := y[2*i]+(1/2)*h*(k[1]+k[2])
end do;
evalf([seq([x[2*i], y[2*i]], i = 0 .. N)])
end proc
## maintenant procedure pour calculer l'erreur (output epsilon(x0,h)
Heunstencil := proc (x0, h)
local g, epsilon;
Heun1((x, y) ->f(x, y),1, h);
Heun2((x, y)->f(x, y), 1, h);
g := (x0, h, i) ->(Heun1( (x, y)->f(x, y) , x0, h)[i][2]-Heun2(proc (x, y)->f(x, y), x0, h)[i][2])^2;
epsilon := (x0, h) ->sqrt(add(g(x0, h, i), i = 1 .. round((1/2)*x0/h)+1)/(round((1/2)*x0/h)+1)) ;
epsilon(x0, h)
end proc
Maintenant si je fais des tests avec x0=1, h=1e-1..1, ca fonctione.
1) Que pensez vous de ce que j'ai fait, car je vais adapter la meme méthode pour Runge Kutta d'ordre 4
2) dans le pdf, l'erreur est defini avec y_i^(2h)-y_(2i)^(h), ce que je fais est-il juste ( je pense que c'est bon) mais toute est relative dans la vie
3) comment peut on tracer en log log la variation de epsilon(x0,h) en fonction de h, pour un x0 fixe.
4) Merci de me donner une idée sur la question Numero 1 ( apropos de l'ordre de p) de l'erreur.
Merci d'avance, et heureux d'echanger des idées avec vous.
-----