heigh_ode.m

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

function [x,y]=heigh_ode(ode_fun, x0, y0, h, N)
% 求解高阶微分方程，其中
%    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; A=[zeros(m-1,1) eye(m-1)];
for n=1:N
    x(n+1)=x(n)+h;
    k1=h*[A*y(:,n); feval(ode_fun, x(n),y(:,n))];
    k2=h*[A*(y(:,n)+1/2*k1); feval(ode_fun, x(n)+1/2*h,y(:,n)+1/2*k1)];
    k3=h*[A*(y(:,n)+1/2*k2); feval(ode_fun, x(n)+1/2*h,y(:,n)+1/2*k2)];
    k4=h*[A*(y(:,n)+k3); feval(ode_fun, x(n)+h,y(:,n)+k3)];
    y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4);
end
y=y(1,:);