Vierte Ordnung Runge–Kutta-Verfahren (RK4) klappt nach ein paar Iterationen
Ich versuche zu lösen:
x' = 60*x - 0.2*x*y;
y' = 0.01*x*y - 100* y;
über die vierte Ordnung Runge-Kutta-Algorithmus.
Ausgangspunkte: x(0) = 8000, y(0) = 300
Bereich: [0,15]
Hier die komplette Funktion:
function [xx yy time r] = rk4_m(x,y,step)
A = 0;
B = 15;
h = step;
iteration=0;
t = tic;
xh2 = x;
yh2 = y;
rr = zeros(floor(15/step)-1,1);
xx = zeros(floor(15/step)-1,1);
yy = zeros(floor(15/step)-1,1);
AA = zeros(1, floor(15/step)-1);
while( A < B)
A = A+h;
iteration = iteration + 1;
xx(iteration) = x;
yy(iteration) = y;
AA(iteration) = A;
[x y] = rkstep(x,y,h);
for h2=0:1
[xh2 yh2] = rkstep(xh2,yh2,h/2);
end
r(iteration)=abs(y-yh2);
end
time = toc(t);
xlabel('Range');
ylabel('Value');
hold on
plot(AA,xx,'b');
plot(AA,yy,'g');
plot(AA,r,'r');
fprintf('Solution:\n');
fprintf('x: %f\n', x);
fprintf('y: %f\n', y);
fprintf('A: %f\n', A);
fprintf('Time: %f\n', time);
end
function [xnext, ynext] = rkstep(xcur, ycur, h)
kx1 = f_prim_x(xcur,ycur);
ky1 = f_prim_y(xcur,ycur);
kx2 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx1);
kx3 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx2);
kx4 = f_prim_x(xcur+h,ycur+h*kx3);
ky2 = f_prim_y(xcur+0.5*h*ky1,ycur+0.5*h);
ky3 = f_prim_y(xcur+0.5*h*ky2,ycur+0.5*h);
ky4 = f_prim_y(xcur+h*ky2,ycur+h);
xnext = xcur + (1/6)*h*(kx1 + 2*kx2 + 2*kx3 + kx4);
ynext = ycur + (1/6)*h*(ky1 + 2*ky2 + 2*ky3 + ky4);
end
function [fx] = f_prim_x(x,y)
fx = 60*x - 0.2*x*y;
end
function [fy] = f_prim_y(x,y)
fy = 0.01*x*y - 100*y;
end
Und ich bin es durch ausführen von: [xx yy time] = rk4_m(8000,300,10)
Das problem ist, dass alles stürzt nach 2-3 Iterationen Rückkehr nutzlose Ergebnisse. Was mache ich falsch? Oder ist diese Methode nicht geeignet für diese Art Gleichung?
Die Semikolons sind absichtlich weggelassen.
Sieht aus wie ich nicht zahlen, die Aufmerksamkeit auf die tatsächlichen h
Größe. Es funktioniert jetzt! Danke!
Du musst angemeldet sein, um einen Kommentar abzugeben.
Sieht aus wie eine form des Lotka-Volterra-Gleichung?
Ich bin mir nicht sicher, ob wenn Ihr ursprünglicher Zustand ist
[300;8000]
oder[8000;300]
(geben Sie es in beide Richtungen oben), aber egal, Sie haben ein schwingungsfähiges system, dass Sie versuchen, zu integrieren, mit einem großen festen Zeit Schritt, die (viel) größer als die Periode der Schwingung. Dies ist der Grund, warum Sie Ihre Fehler explodiert. Wenn Sie versuchen, die Steigerungn
(sagen wir1e6
), werden Sie feststellen, dass schließlich erhalten Sie eine stabile Lösung (vorausgesetzt, dass der Runge-Kutta-Implementierung ist, sonst korrigieren).Gibt es einen Grund, warum Sie nicht mit Hilfe von Matlab-s-builtin ODE-Solver, z.B.
ode45
oderode15s
?Werden Sie feststellen, dass adaptive Schrittweite Löser sind viel effizienter für diese Arten von oszillatorischen Problemen. Denn Ihr system hat eine so hohe Frequenz und scheint ziemlich steif, schlage ich vor, dass Sie auch einen Blick auf das, was
ode15s
gibt und/oder passen Sie die'AbsTol'
und'RelTol'
Optionen mitodeset
.[8000,300]
es war ein Tippfehler in der Frage