rk_four_step.m

微分方程数值解法源代码

function data = RK_four_step(x0,y0,h,N)
%函数功能：四阶龙格－库塔法
%输入参数：（x0，y0）－－初值；h－－步长；N－－计算次数.
%返 回 值：解曲线上的点，格式为二维列向量，第一列为自变量，第二列为函数值.

syms x y;
diff_f = input('请输入微分函数dy/dx:');

data(1,1) = x0;
data(1,2) = y0;

for i = 1:N
    
    %计算K1
    x = x0;
    y = y0;
    K1 = eval(diff_f);
    
    %计算K2
    x = x0 + 0.5 * h;
    y = y0 + 0.5 * h * K1;
    K2 = eval(diff_f);
    
    %计算K3
    x = x0 + 0.5 * h;
    y = y0 + 0.5 * h * K2; 
    K3 = eval(diff_f);
    
    %计算K4
    x = x0 + h;
    y = y0 + h * K3;
    K4 = eval(diff_f);
    
    %计算y1
    y = y0 + h * (K1 + 2 * K2 + 2 * K3 + K4) / 6;
    
    %保存数据
    data(i+1,1) = x;
    data(i+1,2) = y;
    
    %更新数据
    x0 = x;
    y0 = y;
end