Лецкия 9: Метод на Ойлер за решаване на обикновени диференциални уравнения

Задача1:
y’ = 1 + (y/x), x ∈[1,2]
y(1) = 0
(Subscript[y, точно])(x)=

Проверяваме точното решение в края на интервала (там имаме натрупана най-голяма грешка)

2*Log[2.]

f[x_,y_]:=1+(y/x)
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);
For[x=x0; y=y0; i=0, i≤n,i++,Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y]];
y=y+h*f[x,y];
x=x+h]

Намесваме и точното решение:

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);
For[x=x0; y=y0; i=0, i≤n,i++,Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];
y=y+h*f[x,y];
x=x+h]

Увеличаваме броя на точките n = 40:

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=40; h=((xend-x0)/n);
For[x=x0; y=y0; i=0, i≤n,i++,Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];
y=y+h*f[x,y];
x=x+h]

Визуализираме решението от точки и точното решение. За тази цел е нужно да въведем масив от точки.

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]]
TableForm[xy]
ListPlot[xy,Joined->True]
Plot[yt[t],{t,x0,xend}]
Show[%,%%]

Малко украшения:

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%,%%]

Намаляваме стъпката (увеличаваме броя на точките)

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=40; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%,%%]

Задача:
y’ = y - 2 sin x, x∈ [0,0.5]
y(0) = 1
n = 5
(Subscript[y, т]) =sin x + cos x

f[x_,y_]:=y - 2* Sin[x]; yt[x_]:=Sin[x] + Cos[x];
x0=0.; y0 = 1; xend =0.5; n=5; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%,%%]

Увеличаваме интервала, в който интегрираме:

f[x_,y_]:=y - 2* Sin[x]; yt[x_]:=Sin[x] + Cos[x];
x0=0.; y0 = 1; xend =Pi; n=5; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%%,%]

Добавяме повече точки и виждаме как и графиките ни се приближават:

f[x_,y_]:=y - 2* Sin[x]; yt[x_]:=Sin[x] + Cos[x];
x0=0.; y0 = 1; xend =Pi; n=50; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%%,%]

f[x_,y_]:=y - 2* Sin[x]; yt[x_]:=Sin[x] + Cos[x];
x0=0.; y0 = 1; xend =Pi; n=500; h=((xend-x0)/n);
xy={{x0,y0}};
For[x=x0; y=y0; i=0, i≤n,i++,(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y], " yt=", yt[x]];*)
y=y+h*f[x,y];
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy];
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%%,%]

Задача за система от 2 ОДУ

y=.;x=.;z=.
DSolve[{y'[x]==x-y[x] z[x], z'[x] == 5x -z[x], y[0]==1, z[0]==z0},{y[x],z[x]},{x,0,0.4}]

f[x_,y_,z_]:=x-y*z
g[x_,y_,z_]:= 5x-z
y0=1;z0=0;x0=0.; xend=0.4; n=4; h = ((xend-x0)/n);
For[x=x0;y=y0;z=z0; i=0,i≤n,i++,Print[ " i=", i," x=", x," y=", y," z=", z," f(x,y,z)=", f[x,y,z]," g(x,y,z)=", g[x,y,z]];
y1=y+h*f[x,y,z];
z=z+h*g[x,y,z];
y=y1;
x=x+h]

Диференциална задача от втори ред - преминаваме чрез полагане към решаване на система:

f[x_,y_,z_]:=z
g[x_,y_,z_]:= 2z + x y + x Sin[x]
y0=1;z0=-1;x0=0.; xend=1; n=5; h = ((xend-x0)/n);
For[x=x0;y=y0;z=z0; i=0,i≤n,i++,Print[ " i=", i," x=", x," y=", y," z=", z," f(x,y,z)=", f[x,y,z]," g(x,y,z)=", g[x,y,z]];
y1=y+h*f[x,y,z];
z=z+h*g[x,y,z];
y=y1;
x=x+h]

Опитваме се да решим задачата точно с помощта на вградените функции в Mathematica:

DSolve[{y''[x]==2 y'[x]+x*y[x]+x*Sin[x],y[0]==1,y'[0]==-1},y[x],x]

Не можем да го сметнем точно, затова намаляваме стъпката и смятаме числено с по-голяма точност:

f[x_,y_,z_]:=z
g[x_,y_,z_]:= 2z + x y + x Sin[x]
y0=1;z0=-1;x0=0.; xend=1; n=50; h = ((xend-x0)/n);
For[x=x0;y=y0;z=z0; i=0,i≤n,i++,Print[ " i=", i," x=", x," y=", y," z=", z," f(x,y,z)=", f[x,y,z]," g(x,y,z)=", g[x,y,z]];
y1=y+h*f[x,y,z];
z=z+h*g[x,y,z];
y=y1;
x=x+h]