龍格庫塔求解常微分方程

fun.m添加微分方程，經過RK遞推下一時刻的函數值

 

主函數以下：

n=90;

x=zeros(1,n+1);
y=zeros(1,n+1);
x(1)=0;
y(1) =1;   %初值
h=0.1;

for i =1:n
    
    x(i+1) = x(i) + h;
    
    y(i+1) = RK(@fun,x(i),y(i),h);
    
    
end

plot(x,y,'-o')

 

fun.m以下：y' = 2*(1-y/20)*y -x;  

 

function dy= fun(x,y)

dy  = 2*(1-y/20)*y -x;   

　

RK.m以下：

function y = RK(F_xy,x,y,h)


k_1 = F_xy(x,y);

k_2 = F_xy(x+0.5*h,y+0.5*h*k_1);

k_3 = F_xy((x+0.5*h),(y+0.5*h*k_2));

k_4 = F_xy((x+h),(y+k_3*h));


y = y + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  

　　

 

 結果以下：