Лецкия 10: Метод на Рунге-Кута за решаване на обикновени диференциални уравнения

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

Реализираме метод на Рунге-Кута с то4ност O((h^3))

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

Оценка на то4ността:

Локална:

(h^3)

Глобална (в край на интервала)

(h^2)

Правим само две стъпки:

f[x_,y_]:=1+(y/x); yt[x_]:=x*Log[x];
x0=1.; y0 = 0; xend =2; n=2; h=((xend-x0)/n);param=0.5;
xy={{x0,y0}};
p1=param; p2=1-param; α=(1/(2 param));β=α;
Print["  RK O((h^3))  ", " p1 = ", p1, " p2 = ", p2, " α = ", α, " β = ",β]
For[x=x0; y=y0; i=0, i≤n,i++,
k1 = h*f[x,y];
k2=h*f[x+α h, y+ β*k1];Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y]," k1 = ",k1," k2 = ",k2,  " yt=", yt[x]];
y=y+p1*k1+p2*k2;
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=10; h=((xend-x0)/n);param=0.5;
xy={{x0,y0}};
p1=param; p2=1-param; α=(1/(2 param));β=α;
Print["  RK O((h^3))  ", " p1 = ", p1, " p2 = ", p2, " α = ", α, " β = ",β]
For[x=x0; y=y0; i=0, i≤n,i++,
k1 = h*f[x,y];
k2=h*f[x+α h, y+ β*k1];Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y]," k1 = ",k1," k2 = ",k2,  " yt=", yt[x]];
y=y+p1*k1+p2*k2;
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_]:=(Sqrt[x])+(y/x); yt[x_]:=2 (-x+(x^(3/2)));
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);param=0.5;
xy={{x0,y0}};
p1=param; p2=1-param; α=(1/(2 p2));β=α;
Print["  RK O((h^3))  ", " p1 = ", p1, " p2 = ", p2, " α = ", α, " β = ",β]
For[x=x0; y=y0; i=0, i≤n,i++,
k1 = h*f[x,y];
k2=h*f[x+α h, y+ β*k1];(*Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y]," k1 = ",k1," k2 = ",k2,  " yt=", yt[x]];*)
y=y+p1*k1+p2*k2;
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy]
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%,%%]

Пробваме да решим точно:

x=.
y=.
DSolve[{y'[x]==(Sqrt[x])+((y[x])/x),y[1]==0},y[x],x ]

x=.
y=.
s=NDSolve[{y'[x]==(Sqrt[x])+((y[x])/x),y[1]==0},y[x],{x,1,2} ]

Plot[Evaluate[y[x]/.s],{x,1,2},PlotRange->All]

Метод Рунге-Кута с точност O((h^5))

f[x_,y_]:=(Sqrt[x])+(y/x); yt[x_]:=2 (-x+(x^(3/2)));
x0=1.; y0 = 0; xend =2; n=4; h=((xend-x0)/n);
xy={{x0,y0}};
p1=param; p2=1-param; α=(1/(2 p2));β=α;
Print["  RK O((h^3))  "]
For[x=x0; y=y0; i=0, i≤n,i++,
k1 = h*f[x,y];
k2=h*f[x+ (h/2), y+(k1/2)];
k3=h*f[x+ (h/2), y+(k2/2)];
k4=h*f[x+ h, y+k3];Print[" i=", i, " x=", x, " y=", y, " f=", f[x,y]," k1 = ",k1," k2 = ",k2," k3 = ",k3," k4 = ",k4,  " yt=", yt[x]];
y=y+(1/6)(k1+2*k2+2 k3+k4);
x=x+h;
AppendTo[xy,{x,y}]];
TableForm[xy]
ListPlot[xy,Joined->True, PlotStyle->Orange];
Plot[yt[t],{t,x0,xend}];
Show[%,%%]