Решаване на системи ОДУ по метода на Ойлер

задача 4 в). стр. 36 от Ръководство по ЧМ2

ВНИМАНИЕ: ИЗПОЛВАМЕ ВЕЧЕ ПРОМЕНЕНАТА СТОЙНОСТ НА Y

x0=0.; xend=0.4;n=4;
y0=1;z0=0;
h=((xend-x0)/n);
f[x_,y_,z_]:=x-y z;
g[x_,y_,z_]:=5x-z;
(*yt[x_]:=x Log[x];  въвеждаме точното решение *)
x=x0;
y=y0;
z=z0;
(*points = {{x,y}};*)
Print["i = ",0, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
y=y+h*f[x,y,z];
z=z+h*g[x,y,z];
x=x+h;
(* ВНИМАНИЕ: ИЗПОЛВАМЕ ВЕЧЕ ПРОМЕНЕНАТА СТОЙНОСТ НА Y *)
(*AppendTo[points,{x,y}];*)
Print["i = ",i, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
]
(*grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.01]}];
Show[grfyt,grfp]*)

запазваме старите стойности на y и z в отделни променливи:

x0=0.; xend=0.4;n=4;
y0=1;z0=0;
h=((xend-x0)/n);
f[x_,y_,z_]:=x-y z;
g[x_,y_,z_]:=5x-z;
(*yt[x_]:=x Log[x];  въвеждаме точното решение *)
x=x0;
y=y0;
z=z0;
(*points = {{x,y}};*)
Print["i = ",0, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
y1=y+h*f[x,y,z];
z1=z+h*g[x,y,z];
x=x+h;
y=y1; z=z1;
(*AppendTo[points,{x,y}];*)
Print["i = ",i, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
]
(*grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];
grfp=ListPlot[points, PlotStyle->{Black, PointSize[0.01]}];
Show[grfyt,grfp]*)

задача 4 г). стр. 36 от Ръководство по ЧМ2 - самостоятелна работа

задача 4 в). стр. 36 от Ръководство по ЧМ2 - ВИЗУАЛИЗАЦИЯ

x0=0.; xend=0.4;n=4;
y0=1;z0=0;
h=((xend-x0)/n);
f[x_,y_,z_]:=x-y z;
g[x_,y_,z_]:=5x-z;
(*yt[x_]:=x Log[x];  въвеждаме точното решение *)
x=x0;
y=y0;
z=z0;
pointsy = {{x,y}};
pointsz = {{x,z}};
Print["i = ",0, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
y1=y+h*f[x,y,z];
z1=z+h*g[x,y,z];
x=x+h;
y=y1; z=z1;
AppendTo[pointsy,{x,y}];
AppendTo[pointsz,{x,z}];
Print["i = ",i, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z] ]
]
(*grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Red];*)
grfpy=ListPlot[pointsy, PlotStyle->{Blue, PointSize[0.03]}]
grfpz=ListPlot[pointsz, PlotStyle->{Red, PointSize[0.03]}]
Show[grfpy,grfpz, PlotRange->All]