几次迭代后,四阶Runge–Kutta方法(RK4)崩溃 [英] Fourth-order Runge–Kutta method (RK4) collapses after a few iterations
问题描述
我正在尝试解决:
x' = 60*x - 0.2*x*y;
y' = 0.01*x*y - 100* y;
使用四阶Runge-Kutta算法.
using the fourth-order Runge-Kutta algorithm.
起点: x(0)= 8000,y(0)= 300
范围: [0,15]
这是完整的功能:
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
我通过执行以下命令来运行它: [xx yy time] = rk4_m(8000,300,10)
And I'm running it by executing: [xx yy time] = rk4_m(8000,300,10)
问题在于,经过2-3次迭代后,所有内容都会崩溃,并返回无用的结果.我究竟做错了什么?还是这种方法不适用于这种方程式?
The problem is that everything collapses after 2-3 iterations returning useless results. What am I doing wrong? Or is just this method not appropriate for this kind equation?
有意省略分号.
好像我没有注意实际的 h
大小.现在可以使用了!谢谢!
Looks like I didn't pay attention to actual h
size. It works now! Thanks!
推荐答案
看起来像某种形式的Lotka-Volterra方程?
Looks like some form of the Lotka-Volterra equation?
我不确定您的初始条件是 [300; 8000]
还是 [8000; 300]
(您在上述两种方式中都指定了),但是无论,您有一个振荡系统,您正试图与一个较大的固定时间步长(该步长比振荡周期大得多)相集成.这就是为什么您的错误会爆炸的原因.如果您尝试增加 n
(例如, 1e6
),则会发现最终您将获得一个稳定的解决方案(假设您的Runge-Kutta实现正确)
I'm not sure if if your initial condition is [300;8000]
or [8000;300]
(you specify it both ways above), but regardless, you have an oscillatory system that you're trying to integrate with a large fixed time step that is (much) greater than the period of oscillation. This is why your error explodes. If you try increasing n
(say, 1e6
), you'll find that eventually you'll get a stable solution (assuming that your Runge-Kutta implementation is otherwise correct).
为什么没有使用Matlab的内置ODE求解器,例如 ode45
或
Is there a reason why you're not using Matlab's builtin ODE solvers, e.g. ode45
or ode15s
?
function ode45demo
[t,y]=odeode45(@f,[0 15],[300;8000]);
figure;
plot(t,y);
function ydot=f(t,y)
ydot(1,1) = 60*y(1) - 0.2*y(1)*y(2);
ydot(2,1) = 0.01*y(1)*y(2) - 100*y(2);
您会发现,自适应步长求解器对于这类振荡问题更为有效.因为您的系统频率很高,而且看起来比较僵硬,所以我建议您还查看 ode15s
给出的内容和/或调整'AbsTol'
和'带有
odeset
选项>.
You'll find that adaptive step size solvers are much more efficient for these types of oscillatory problems. Because your system has such a high frequency and seems rather stiff, I suggest that you also look at what ode15s
gives and/or adjust the 'AbsTol'
and 'RelTol'
options with odeset
.
这篇关于几次迭代后,四阶Runge–Kutta方法(RK4)崩溃的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!