Модифициран Ойлер за системи ОДУ

съставяме задача с точно решение:

x=.
yt[x_]:= x+Sin[x]
zt[x_]:= 2x - Cos[x]
yt'[x]
zt'[x]

yt[0]

zt[0]

Задачата е:
y’ = -z + 2x + 1
z’ = y + z + 2 - 3x + cos x
y(0) = 0
z(0) = -1
x ∈ [0, 1]
точно решение: y*(x) = x + sin x, z*(x) = 2x - cos x

самостоятелна проверка на точното решение и проба за получаване с DSolve

програмен код

метод на Ойлер

x0=0.; xend=1;n=5;
y0=0;z0=-1;
h=((xend-x0)/n);
f[x_,y_,z_]:=-z +2x +1;
g[x_,y_,z_]:=y + z +2 -3x + Cos[x];
yt[x_]:= x+Sin[x]
zt[x_]:= 2x - Cos[x]  (*въвеждаме точното решение *)
x=x0;
y=y0;
z=z0;
pointsy = {{x,y}};
pointsz = {{x,z}};
epsy = Abs[y-yt[x]];
epsz = Abs[z-zt[x]];
Print["i = ",0, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z]," (SuperStar[y]) = ",yt[x], " (SuperStar[z]) = ",zt[x], " (Subscript[ϵ, y]) = ",epsz, " (Subscript[ϵ, z]) = ",epsz]
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}];
epsy = Abs[y-yt[x]];
epsz = Abs[z-zt[x]];
Print["i = ",i, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z]," (SuperStar[y]) = ",yt[x], " (SuperStar[z]) = ",zt[x], " (Subscript[ϵ, y]) = ",epsz, " (Subscript[ϵ, z]) = ",epsz ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Black];
grfzt=Plot[zt[x],{x,x0,xend}, PlotStyle->Green];
grfpy=ListPlot[pointsy, PlotStyle->{Blue, PointSize[0.03]}];
grfpz=ListPlot[pointsz, PlotStyle->{Red, PointSize[0.03]}];
Show[grfpy,grfpz,grfyt,grfzt, PlotRange->All]

модифициран метод на Ойлер

x0=0.; xend=1;n=5;
y0=0;z0=-1;
h=((xend-x0)/n);
f[x_,y_,z_]:=-z +2x +1;
g[x_,y_,z_]:=y + z +2 -3x + Cos[x];
yt[x_]:= x+Sin[x]
zt[x_]:= 2x - Cos[x]  (*въвеждаме точното решение *)
x=x0;
y=y0;
z=z0;
pointsy = {{x,y}};
pointsz = {{x,z}};
epsy = Abs[y-yt[x]];
epsz = Abs[z-zt[x]];
Print["i = ",0, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z]," (SuperStar[y]) = ",yt[x], " (SuperStar[z]) = ",zt[x], " (Subscript[ϵ, y]) = ",epsz, " (Subscript[ϵ, z]) = ",epsz]
For[i=1,i≤ n,i++,
(* тяло на цикъла *)
x12 =x+(h/2) ;y12 = y+(h/2)f[x,y,z];z12 =z+(h/2)*g[x,y,z] ;
y1=y+h*f[x12,y12,z12];
z1=z+h*g[x12,y12,z12];
x=x+h;
y=y1; z=z1;
AppendTo[pointsy,{x,y}];
AppendTo[pointsz,{x,z}];
epsy = Abs[y-yt[x]];
epsz = Abs[z-zt[x]];
Print["i = ",i, " x = ",x, " y = ",y, " z = ",z, " f(x,y) = ",f[x,y,z]," g(x,y) = ",g[x,y,z]," (SuperStar[y]) = ",yt[x], " (SuperStar[z]) = ",zt[x], " (Subscript[ϵ, y]) = ",epsz, " (Subscript[ϵ, z]) = ",epsz ]
]
grfyt=Plot[yt[x],{x,x0,xend}, PlotStyle->Black];
grfzt=Plot[zt[x],{x,x0,xend}, PlotStyle->Green];
grfpy=ListPlot[pointsy, PlotStyle->{Blue, PointSize[0.03]}];
grfpz=ListPlot[pointsz, PlotStyle->{Red, PointSize[0.03]}];
Show[grfpy,grfpz,grfyt,grfzt, PlotRange->All]