r-k4.m

数值分析常用几种方法

function R=r-k4(f,a,b,f0,h)
%f为函数f(x,y);
%a,b为取值范围的左右端点；
%f0为函数初值y(0)，h为步长,n为迭代步数；
n=(b-a)/h;
T=zeros(1,n+1);
Y=zeros(1,n+1);
T=a:h:b;
Y(1)=f0;
for j=1:n
    K1=h*feval('f',T(j),Y(j));
    K2=h*feval('f',T(j)+h/2,Y(j)+K1/2);
    K3=h*feval('f',T(j)+h/2,Y(j)+K2/2);
    K4=h*feval('f',T(j)+h,Y(j)+K3);
    Y(j+1)=Y(j)+(K1+2*K2+2*K3+K4)/6;
end
t=[0:.1:5];
y=1./exp(10*t);
figure(1)
plot(T,Y,'*',T,Y,'r-',t,y,'k-');
hold on
R=[T' Y']