rk.m

龙格库塔程序源代码

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

syms t x x1 x2;
diff_f = input('请输入微分函数dy/dx(列向量):');

dx1 = diff_f(1,1);
dx2 = diff_f(2,1);

data(1,1) = t0(1,1);
data(1,2) = x0(1,1);
data(1,3) = t0(2,1);
data(1,4) = x0(2,1);

for i = 1:N
    t0 = data(i,1);
    x0 = data(i,2);
    x1 = x;
    x2 = data(i,3);
    dx = eval(dx1);
    data(i+1,(1:2)) = RK_sigle_step(t0,x0,h,dx);

    t0 = data(i,3);
    x0 = data(i,4);
    x1 = data(i,2);
    x2 = x;
    dx = eval(dx2);
    data(i+1,(3:4)) = RK_sigle_step(t0,x0,h,dx);
end
plot(data(:,1),data(:,2),'r',data(:,3),data(:,4),'b');