實現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

回答

相關問題