www.pudn.com > MATLABpro.rar > RK.m

%四阶RK法求解常微分方程 
function y=RK(a,b,N,af); 
h=(b-a)/N; 
x(1)=a; 
y(1)=af; 
jqj(1)=af; 
for i=2:N+1 
    K1=f(x(i-1),y(i-1)); 
    K2=f(x(i-1)+h/2,y(i-1)+h*K1/2); 
    K3=f(x(i-1)+h/2,y(i-1)+h*K2/2); 
    K4=f(x(i-1)+h,y(i-1)+h*K3); 
    y(i)=y(i-1)+(K1+2*K2+2*K3+K4)/6; 
    x(i)=x(i-1)+(i-1)*h; 
    jqj(i)=x(i)+exp((-x(i))); 
end 
[x',y',jqj'] 
er=norm(y-jqj,2)/norm(y) 
plot(x',y','r',x',jqj','g'); 
legend('RK法','精确解');