龙格库塔求解常微分方程

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;

结果如下：

原文地址：https://www.cnblogs.com/skylover/p/10160892.html