euler_ts.m

来自「插值matlab源程序,附使用说明,包括拉各朗日插值,牛顿插值,及示例」· M 代码 · 共 32 行

function [x,y]=euler_ts(ydot_fun, x0, y0, h, N)
% 梯形公式，其中
%    ydot_fun ---  一阶微分方程的函数
%    x0, y0   ---  初始条件
%    h        ---  区间步长
%    N        ---  区间的个数
%    x        ---  Xn 构成的向量
%    y        ---  Yn 构成的向量

m=length(y0);x=zeros(1,N+1); y=zeros(m,N+1); 
x(1)=x0; y(:,1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    % 用迭代法求y(n+1)
    z0=y(:,n)+h*feval(ydot_fun, x(n), y(:,n));
    for k=1:10
        z1=y(:,n)+h/2*(feval(ydot_fun, x(n), y(:,n)) ... 
               +feval(ydot_fun, x(n+1), z0));
        z2=y(:,n)+h/2*(feval(ydot_fun, x(n), y(:,n)) ... 
               +feval(ydot_fun, x(n+1), z1));
        if norm(z2-z1)<1e-5
            break;
        end
        z1=z0-((z1-z0).^2)./(z2-2*z1+z0);
        if norm(z1-z0)<1e-5
            break;
        end
        z0=z1;
    end
    y(:,n+1)=z1;
end