实现Runge-Kutta 4

Question

我正在尝试实现Runge 4来求解一个微分方程，并需要一些洞察力。

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

我重复一遍:我不是在找解决方案，我在找人帮我指出我的代码出了什么问题，因为我已经盯着这个看了多久了.

代码的工作原理:在函数声明和for循环之间设置初始值(同样，p和q在问题的这一部分中无关紧要)，for循环在数值上迭代和求解函数。我使用runge 4，如维基百科文章所述。

我编写的代码如下，N=400 (400个步骤)，T=42 (终端时间4pi2)，(x，y)=(1,0) (初始条件)，gamma =1(γ是方程中的一个参数)，(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步骤。最重要的是，您对错误的定义是错误的。您想要检查半径是否为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