实现Runge-Kutta 4

我想实现Runge-Kutta 4来求解一个微分方程并需要一些见解。实现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（γ是方程中的一个参数），（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

下面是功能的MATLAB代码我解决（P1PC）：

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

来源

2015-01-20 malxmusician212

这是一个二维问题？ – patrik 2015-01-20 22:30:19

是的，kx和ky是我对dx/dt，dy/dt的简写，尽管我认为写dx和dy会更有意义lol – malxmusician212 2015-01-20 22:45:53

你把它和'ode45（）'的结果比较了吗？ – ja72 2015-01-20 22:54:48

我看到两个问题。一个是只有在采取N步骤后才能达到结束时间，并且您的代码需要执行N-1步骤。最重要的是，你对错误的定义是错误的。要检查是否半径是一个，因此error=sqrt(x^2+y^2)-1

请参见下面的代码，我用（有点简化）来测试算法

P1PC [email protected](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

来源

2015-01-20 23:25:49 ja72

你真是个好人！非常感谢！ – malxmusician212 2015-01-21 01:14:39

另一个快乐的客户！ – ja72 2015-01-21 03:29:38

实现Runge-Kutta 4

回答

相关问题