kutta.m

来自「用四阶古典RK方法解初值问题： 取h=1/8。每隔8步打印出数值解与真解」· M 代码 · 共 32 行

%四阶龙格库塔方法
function[T,U]=kutta(a,b,h,c)
%计算点的个数
n=round((b-a)/h)+1;
%定义T,U的维数
t=zeros(n,1);
u=zeros(n,1);
%微分方程的初始条件
t(1)=a;
u(1)=c;
%显示微分方程的初始条件
str=sprintf('t0=%g, u0=%g\n',t(1),u(1));
disp(str);
%对每个点用四阶龙格库塔方法计算
for i=1:n-1
    %系数
    k1=f(t(i),u(i));
    k2=f(t(i)+0.5*h,u(i)+0.5*h*k1);
    k3=f(t(i)+0.5*h,u(i)+0.5*h*k2);
    k4=f(t(i)+h,u(i)+h*k3);
    %龙格库塔公式
    u(i+1)=u(i)+h*(k1+2*k2+2*k3+k4)/6;
    t(i+1)=t(i)+h;
    %显示计算结果与误差
    str=sprintf('%d t=%g u=%g orig=%g e=%g',i,t(i+1),u(i+1),orig(t(i+1)),abs(u(i+1)-orig(t(i+1))));
     if rem(i,8)==0;
         disp(str);
     end
 end
 %返回计算结果
 T=t;
 U=u;