rk3.m

来自「常微分方程的龙格库塔法系列程序」· M 代码 · 共 25 行

function rk3(f,u0,a,b,h)
%Runge_Kutta法,s=3
%调用格式:rk3(function_handle,u0,a,b,h)
%test: du/dt=u/2;u0=1;h=0.1;a=0;b=1;
%调用:先定义微分方程f=@(t,u) u/2;rk3(f,1,0,1,0.1)

%2007-03-15
tic;
len=1+(b-a)/h;
u=zeros(1,len);
t=a:h:b;
u(1)=u0;
for i=2:len,
    K1=f(t(i-1),u(i-1));
    K2=f(t(i-1),u(i-1)+h*K1/2);
    K3=f(t(i-1),u(i-1)+3*h*K2/4);
    u(i)=u(i-1)+h*(2*K1+3*K2+4*K3)/9;
end
plot(t,u,'b+');
hold on;
line(t,u,'Color','b');
fun=@(t) exp(t/2);
fplot(fun,[0,1],'r');
error=sqrt(sum(abs(fun(t)-u).^2))
toc