adams_wjl1.m

(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla

function [x,y]=adams_wjl1(odefile,xi,xf,yi,h,varargin)

% odefile is the filename of ode , xi<x <xf,  yi is the initial value, h is
% the step , varargin  is used by odefile

%
%  See also ODE23, ODE45, ODE113, ODE15S, ODE23S,
%           

x = [xi:h:xf];		% Vector of x values
if x(end) ~= xf
   x(end+1) = xf;
end

d = diff(x);
% if you don't have this sentence,change d(i) to h is OK

% Starting values
% tic;

% [a,b] = rk_wjl(odefile,x(1),x(4),h,yi,varargin{:});
[a,b] = rk_wjl(odefile,x(1),x(4),yi,h,varargin{:});
%do not write yi , h, in the wrong direction
y(1:4) = b;
for i = 1:4
   f(i) = feval(odefile,x(i),y(:,i),varargin{:});
end

ynp=0;
ync=0;



% Solution
for i = 4:length(x)-1
   % Predictor
   yn1p = y(i) + d(i)/24 * (55*f(i) - 59*f(i-1) ...
      + 37*f(i-2) - 9*f(i-3));
   yn1pm=yn1p +251/270*(ync-ynp);
   
   %f(i+1) = feval(odefile,x(i+1),y(i+1),varargin{:});
    f(i+1) = feval(odefile,x(i+1),yn1pm,varargin{:});

   % Corrector
   yn1c = y(i) + d(i)/24 * (9*f(i+1) + 19*f(i) ...
      - 5*f(i-1) + f(i-2));
   y(i+1) =yn1c -19/270*(yn1c - yn1p);
   %f(i+1) = feval(odefile,x(i+1),y(i+1),varargin{:});
   ync=yn1c;
   ynp=yn1p;
   
end
% time=toc;
% you use time out the programme