4lkfun.txt

4阶变步长龙格-库塔法解dy=fun(x,y)型方程,根据算法编写。

4阶变步长龙格-库塔法解dy=fun(x,y)型方程
function change_step_RK(fun); 
% 4阶变步长龙格-库塔法解dy=fun(x,y)型方程 
% Example: 
%      fun=inline('cos(x)+sin(y)'); 
%      change_step_RK(fun); 
AbsTol=1e-4; 
h=0.01; 
p=4;    % p是对应的阶数 
x0=0; % x0是起点 
xN=1; % xN是终点 
y0=0; % 初值 
x=x0; 
y=y0; 
n=1; 
p21=2^p-1; 
while x(end)<xN; 
    [x1,y1]=RK_f(fun,x(n),y(n),h); 
    [x2,y2]=RK_f(fun,x(n),y(n),h/2); 
    if abs(y1-y2)/p21<AbsTol;    % 计算步长 
        while abs(y1-y2)/p21<AbsTol; 
            x2=x1; 
            y2=y1; 
            h=2*h; 
            [x1,y1]=RK_f(fun,x(n),y(n),h); 
        end 
    else 
        while abs(y1-y2)/p21>AbsTol; 
            x1=x2; 
            y1=y2; 
            h=h/2; 
            [x2,y2]=RK_f(fun,x(n),y(n),h/2); 
        end 
    end 
    [xa,ya]=RK_f(fun,h,x(n),y(n)); 
    x(n+1)=xa; 
    y(n+1)=ya; 
    n=n+1; 
end 
plot(x,y,'k'); 
function [xa,ya]=RK_f(fun,h,x,y); 
% 4阶龙格-库塔法解下一点的值 
k1=fun(x,y); 
k2=fun(x+h/2,y+h*k1/2); 
k3=fun(x+h/2,y+h*k2/2); 
k4=fun(x+h,y+h*k3); 
xa=x+h; 
ya=y+h*(k1+k2*2+2*k3+k4)/6;