实施 Runge-Kutta 4

Question

我正在尝试实施 Runge-Kutta 4 来求解微分方程，需要一些见解。

我遇到的主要问题是这些值的误差在 0.09 左右，而它应该在 0.0001*10-4 左右，所以我的 rk4 实现有问题，但我不知道我在哪里犯了错误。如果你能指出我的 runge-kutta 实施的方向，那就太好了。（注意，我们能够计算误差是因为解是一个圆，所以我知道终点应该与初始条件（1,0）相同，误差是计算出的端点与(1,0)。此作业用于学习数值方法。

我再说一遍：我不是在寻找解决方案，我是在找人帮我指出我的代码出了什么问题，因为天知道我已经盯着这个看多久了。 ..

代码的工作原理：在函数声明和 for 循环之间设置初始值（同样，p 和 q 在这部分问题中无关紧要），for 循环迭代并以数值方式求解函数。我按照维基百科文章中的描述使用 runge kutta 4。

我写的代码如下，N=400（400步），T=42（4pi2的终端时间），(x,y)=(1,0)（初始条件），gamma = 1（gamma 是等式中的一个参数），并且 (p,q) 与我要问的部分无关。 P1PC是包含微分方程的.m文件，也在下面。

function rk(N,T,x,y,gamma,p,q)
timestep=T/N;
xy=zeros(N,4);
xy(1,:)=[x,y,p,q];
k=zeros(4,2);%k(:,1) is for x, k(:,2) is for y
for index=2:N
    [k(1,1),k(1,2)]=P1PC(gamma,xy(index-1,1),xy(index-1,2));
    [k(2,1),k(2,2)]=P1PC(gamma,xy(index-1,1)+0.5*k(1,1)*timestep,xy(index-1,2)+0.5*k(1,2)*timestep);
    [k(3,1),k(3,2)]=P1PC(gamma,xy(index-1,1)+0.5*k(2,1)*timestep,xy(index-1,2)+0.5*k(2,2)*timestep);
    [k(4,1),k(4,2)]=P1PC(gamma,xy(index-1,1)+k(3,1)*timestep,xy(index-1,2)+k(3,2)*timestep);
    xy(index,1)=xy(index-1,1)+(timestep/6)*(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1));
    xy(index,2)=xy(index-1,2)+(timestep/6)*(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2));
end
plot(xy(:,1),xy(:,2));
error=sqrt((1-xy(N,1))^2+(0-xy(N,2))^2)
xy(N,1)
xy(N,2)
end

这是我正在求解的函数 (P1PC) 的 matlab 代码：

function [kx,ky]=P1PC(gamma,x,y)
    r=x^2+y^2;
    ky=(gamma*x)/(2*pi*r);
    kx=(gamma*(-y))/(2*pi*r);
    end

Answer 1

我看到两个问题。一种是结束时间仅在执行 N 步后到达，而您的代码执行了 N-1 步。最重要的是，你对错误的定义是错误的。您想检查半径是否为一，因此 error=sqrt(x^2+y^2)-1

请参阅下面我用来（稍微简化）来测试算法的代码

P1PC =@(gamma,x,y)[-gamma*y,gamma*x]/(2*pi*(x^2+y^2));
T = 42;
N = 400;
h=T/N;
time=0;
x0=1;
y0=0;
gamma=1;
xy = zeros(N+1,2);
xy(1,:) = [x0,y0];
for i=2:N+1
    k1 = P1PC(gamma, xy(i-1,1),xy(i-1,2));
    k2 = P1PC(gamma, xy(i-1,1)+(h/2)*k1(1), xy(i-1,2)+(h/2)*k1(2));
    k3 = P1PC(gamma, xy(i-1,1)+(h/2)*k2(1), xy(i-1,2)+(h/2)*k2(2));
    k4 = P1PC(gamma, xy(i-1,1)+(h)*k3(1), xy(i-1,2)+(h)*k3(2));
    xy(i,:) = xy(i-1,:) + (h/6)*(k1+2*k2+2*k3+k4);
    time = time + h;
end

plot(xy(:,1),xy(:,2));
disp(['time=',num2str(time)])
error=sqrt(xy(N+1,1)^2+xy(N+1,2)^2)-1;
disp(['error=',num2str(error)])

产生输出

time=42
error=3.0267e-011

实施 Runge-Kutta 4

Implementing Runge-Kutta 4

matlab

numerical-methods