24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

最怜宵舞

木虫 (正式写手)

应助: 1 (幼儿园)
金币: 1711.4
散金: 114
红花: 16
沙发: 6
帖子: 420
在线: 24小时
虫号: 4171587
注册: 2015-10-25
性别: MM
专业: 计算数学与科学工程计算

[求助] ode45已有4人参与

一道常微分求解的题目，要求用四阶龙格库塔法R-K，求Matlab中的ode45算法的内部算法程序代码，或者思路。

发自小木虫Android客户端

回复此楼

» 猜你喜欢

Bioresource Technology期刊，第一次返修的时候被退回好几次了已经有6人回复
2025冷门绝学什么时候出结果已经有4人回复
真诚求助：手里的省社科项目结项要求主持人一篇中文核心，有什么渠道能发核心吗已经有8人回复
寻求一种能扛住强氧化性腐蚀性的容器密封件已经有5人回复
论文投稿，期刊推荐已经有6人回复
请问哪里可以有青B申请的本子可以借鉴一下。已经有4人回复
孩子确诊有中度注意力缺陷已经有14人回复
请问下大家为什么这个铃木偶联几乎不反应呢已经有5人回复
请问有评职称，把科研教学业绩算分排序的高校吗已经有5人回复
天津工业大学郑柳春团队欢迎化学化工、高分子化学或有机合成方向的博士生和硕士生加入已经有4人回复

1楼2016-01-07 15:41:25

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

wurongjun

专家顾问 (职业作家)

专家经验: +831
数学EPI: 9
应助: 791 (博后)
贵宾: 0.308
金币: 24607.5
散金: 310
红花: 75
帖子: 3004
在线: 881.2小时
虫号: 1368482
注册: 2011-08-14
性别: GG
专业: 计算数学与科学工程计算
管辖: 数学

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
最怜宵舞: 金币+10, ★★★很有帮助 2016-01-08 15:07:38

下面是ode45的源代码:
function varargout = ode45(ode,tspan,y0,options,varargin)
%ODE45  Solve non-stiff differential equations, medium order method.
% [T,Y] = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the
% system of differential equations y' = f(t,y) from time T0 to TFINAL with
% initial conditions Y0. Function ODEFUN(T,Y) must return a column vector
% corresponding to f(t,y). Each row in the solution array Y corresponds to
% a time returned in the column vector T. To obtain solutions at specific
% times T0,T1,...,TFINAL (all increasing or all decreasing), use
% TSPAN = [T0 T1 ... TFINAL].
%
% [T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
% integration properties replaced by values in OPTIONS, an argument created
% with the ODESET function. See ODESET for details. Commonly used options
% are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
% of absolute error tolerances 'AbsTol' (all components 1e-6 by default).
%
% [T,Y] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2...) passes the additional
% parameters P1,P2,... to the ODE function as ODEFUN(T,Y,P1,P2...), and to
% all functions specified in OPTIONS. Use OPTIONS = [] as a place holder if
% no options are set.
%
% ODE45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
% nonsingular. Use ODESET to set the 'Mass' property to a function MASS if
% MASS(T,Y) returns the value of the mass matrix. If the mass matrix is
% constant, the matrix can be used as the value of the 'Mass' option. If
% the mass matrix does not depend on the state variable Y and the function
% MASS is to be called with one input argument T, set 'MStateDependence' to
% 'none'. ODE15S and ODE23T can solve problems with singular mass matrices.
%
% [T,Y,TE,YE,IE] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS...) with the 'Events'
% property in OPTIONS set to a function EVENTS, solves as above while also
% finding where functions of (T,Y), called event functions, are zero. For
% each function you specify whether the integration is to terminate at a
% zero and whether the direction of the zero crossing matters. These are
% the three vectors returned by EVENTS: [VALUE,ISTERMINAL,DIRECTION] =
% EVENTS(T,Y). For the I-th event function: VALUE(I) is the value of the
% function, ISTERMINAL(I)=1 if the integration is to terminate at a zero of
% this event function and 0 otherwise. DIRECTION(I)=0 if all zeros are to
% be computed (the default), +1 if only zeros where the event function is
% increasing, and -1 if only zeros where the event function is
% decreasing. Output TE is a column vector of times at which events
% occur. Rows of YE are the corresponding solutions, and indices in vector
% IE specify which event occurred.
%
% SOL = ODE45(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
% used with DEVAL to evaluate the solution at any point between T0 and
% TFINAL. The steps chosen by ODE45 are returned in a row vector SOL.x.
% For each I, the column SOL.y(:,I) contains the solution at SOL.x(I).
% If events were detected, SOL.xe is a row vector of points at which events
% occurred. Columns of SOL.ye are the corresponding solutions, and indices
% in vector SOL.ie specify which event occurred. If a terminal event has
% been detected, SOL.x(end) contains the end of the step at which the event
% occurred. The exact point of the event is reported in SOL.xe(end).
%
% Example
%       [t,y]=ode45(@vdp1,[0 20],[2 0]);
%       plot(t,y(:,1));
%    solves the system y' = vdp1(t,y), using the default relative error
%    tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
%    component, and plots the first component of the solution.
%
% See also
%    other ODE solvers: ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB
%    options handling:    ODESET, ODEGET
%    output functions:    ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT
%    evaluating solution:  DEVAL
%    ODE examples:       RIGIDODE, BALLODE, ORBITODE
%
% NOTE:
%    The interpretation of the first input argument of the ODE solvers and
%    some properties available through ODESET have changed in this version
%    of MATLAB. Although we still support the v5 syntax, any new
%    functionality is available only with the new syntax. To see the v5
%    help, type in the command line
%       more on, type ode45, more off

% NOTE:
%    This portion describes the v5 syntax of ODE45.
%
% [T,Y] = ODE45('F',TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the
% system of differential equations y' = F(t,y) from time T0 to TFINAL with
% initial conditions Y0.  'F' is a string containing the name of an ODE
% file.  Function F(T,Y) must return a column vector.  Each row in
% solution array Y corresponds to a time returned in column vector T.  To
% obtain solutions at specific times T0, T1, ..., TFINAL (all increasing
% or all decreasing), use TSPAN = [T0 T1 ... TFINAL].
%
% [T,Y] = ODE45('F',TSPAN,Y0,OPTIONS) solves as above with default
% integration parameters replaced by values in OPTIONS, an argument
% created with the ODESET function.  See ODESET for details.  Commonly
% used options are scalar relative error tolerance 'RelTol' (1e-3 by
% default) and vector of absolute error tolerances 'AbsTol' (all
% components 1e-6 by default).
%
% [T,Y] = ODE45('F',TSPAN,Y0,OPTIONS,P1,P2,...) passes the additional
% parameters P1,P2,... to the ODE file as F(T,Y,FLAG,P1,P2,...) (see
% ODEFILE).  Use OPTIONS = [] as a place holder if no options are set.
%
% It is possible to specify TSPAN, Y0 and OPTIONS in the ODE file (see
% ODEFILE).  If TSPAN or Y0 is empty, then ODE45 calls the ODE file
% [TSPAN,Y0,OPTIONS] = F([],[],'init') to obtain any values not supplied
% in the ODE45 argument list.  Empty arguments at the end of the call list
% may be omitted, e.g. ODE45('F').
%
% ODE45 can solve problems M(t,y)*y' = F(t,y) with a mass matrix M that is
% nonsingular.  Use ODESET to set Mass to 'M', 'M(t)', or 'M(t,y)' if the
% ODE file is coded so that F(T,Y,'mass') returns a constant,
% time-dependent, or time- and state-dependent mass matrix, respectively.
% The default value of Mass is 'none'. ODE15S and ODE23T can solve problems
% with singular mass matrices.
%
% [T,Y,TE,YE,IE] = ODE45('F',TSPAN,Y0,OPTIONS) with the Events property in
% OPTIONS set to 'on', solves as above while also locating zero crossings
% of an event function defined in the ODE file.  The ODE file must be
% coded so that F(T,Y,'events') returns appropriate information.  See
% ODEFILE for details.  Output TE is a column vector of times at which
% events occur, rows of YE are the corresponding solutions, and indices in
% vector IE specify which event occurred.
%
% See also ODEFILE

% ODE45 is an implementation of the explicit Runge-Kutta (4,5) pair of
% Dormand and Prince called variously RK5(4)7FM, DOPRI5, DP(4,5) and DP54.
% It uses a "free" interpolant of order 4 communicated privately by
% Dormand and Prince.  Local extrapolation is done.

% Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
% M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997.

% Mark W. Reichelt and Lawrence F. Shampine, 6-14-94
% Copyright 1984-2001 The MathWorks, Inc.
% $Revision: 5.70 $  $Date: 2001/02/16 16:19:59 $

true = 1;
false = ~true;

stats = struct('nsteps',0,'nfailed',0,'nfevals',0,...
            'npds',0,'ndecomps',0,'nsolves',0);

% Check inputs
if nargin < 4
  options = [];
  if nargin < 3
y0 = [];
if nargin < 2
   tspan = [];
   if nargin < 1
      error('Not enough input arguments.  See ODE45.');
   end
end
  end
end

FcnHandlesUsed = isa(ode,'function_handle');
soloutRequested = (FcnHandlesUsed & (nargout==1));

% Set solver arguments
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, ...
      options, atol, rtol, threshold, normcontrol, normy, hmax, htry, htspan]  ...
      = odearguments(FcnHandlesUsed,'ode45', ode, tspan, y0, options,  ...
                     soloutRequested, varargin);
stats.nfevals = stats.nfevals + 1;

% Handle the output
if nargout > 0
  outfun = odeget(options,'OutputFcn',[],'fast');
else
  outfun = odeget(options,'OutputFcn',@odeplot,'fast');
end
if isempty(outfun)
  haveoutfun = false;
else
  haveoutfun = true;
  outputs = odeget(options,'OutputSel',1:neq,'fast');
  if isa(outfun,'function_handle')
outputArgs = [{''},varargin];
outputArgs1 = varargin;
  else    % With v5 syntax do not pass additional input arguments to outfun
outputArgs = {};
outputArgs1 = {};
  end
end
if soloutRequested
  refine = 1;
else
  refine = odeget(options,'Refine',4,'fast');
end
if ntspan > 2
  outflag = 1;                               % output only at tspan points
elseif refine <= 1
  outflag = 2;                               % computed points, no refinement
else
  outflag = 3;                               % computed points, with refinement
  S = (1:refine-1)' / refine;
end
printstats = strcmp(odeget(options,'Stats','off','fast'),'on');

% Handle the event function
[haveeventfun,eventfun,eventargs,valt,teout,yeout,ieout] = ...
odeevents(FcnHandlesUsed,ode,t0,y0,options,varargin);

% Handle the mass matrix
[Mtype, Mfun, Margs, M] =  odemass(FcnHandlesUsed,ode,t0,y0,options,varargin);
if Mtype>0  % non-trivial mass matrix
  Msingular = odeget(options,'MassSingular','no','fast');
  if strcmp(Msingular,'maybe')
warning(['This solver assumes MassSingular is ''no''.  See ODE15S or ' ...
         'ODE23T.']);
  elseif strcmp(Msingular,'yes')
   error(['MassSingular cannot be ''yes'' for this solver.  See ODE15S '...
         ' or ODE23T.']);
  end
  if Mtype == 1
[L,U] = lu(M);
  else
L = [];
U = [];
  end
  [ode,args] = odemassexplicit(FcnHandlesUsed,Mtype,ode,args,Mfun,Margs,L,U);

  % Recompute f0, this time fun takes care about the mass matrix
  f0 = feval(ode,t0,y0,args{:});
  stats.nfevals = stats.nfevals + 1;
end

t = t0;
y = y0;

% Allocate memory if we're generating output.
if nargout > 0
  if ntspan > 2                      % output only at tspan points
tout = zeros(ntspan,1);
yout = zeros(ntspan,neq);
  else                               % alloc in chunks
chunk = min(max(100,50*refine),floor((2^13)/neq));
tout = zeros(chunk,1);
yout = zeros(chunk,neq);
if soloutRequested
   f3d = zeros(chunk,neq,7);
end
  end
  nout = 1;
  tout(nout) = t;
  yout(nout,

= y.';
end

% Initialize method parameters.
pow = 1/5;
A = [1/5; 3/10; 4/5; 8/9; 1; 1];
B = [
1/5       3/40 44/45 19372/6561    9017/3168    35/384
0          9/40 -56/15  -25360/2187    -355/33       0
0          0    32/9 64448/6561    46732/5247    500/1113
0          0    0    -212/729       49/176       125/192
0          0    0    0             -5103/18656    -2187/6784
0          0    0    0             0             11/84
0          0    0    0             0             0
];
E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];
f = zeros(neq,7);
hmin = 16*eps*abs(t);
if isempty(htry)
  % Compute an initial step size h using y'(t).
  absh = min(hmax, htspan);
  if normcontrol
rh = (norm(f0) / max(normy,threshold)) / (0.8 * rtol^pow);
  else
rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
  end
  if absh * rh > 1
absh = 1 / rh;
  end
  absh = max(absh, hmin);
else
  absh = min(hmax, max(hmin, htry));
end
f(:,1) = f0;

% Initialize the output function.
if haveoutfun
  feval(outfun,[t tfinal],y(outputs),'init',outputArgs1{:});
end

% THE MAIN LOOP

done = false;
while ~done

  % By default, hmin is a small number such that t+hmin is only slightly
  % different than t.  It might be 0 if t is 0.
  hmin = 16*eps*abs(t);
  absh = min(hmax, max(hmin, absh)); % couldn't limit absh until new hmin
  h = tdir * absh;

  % Stretch the step if within 10% of tfinal-t.
  if 1.1*absh >= abs(tfinal - t)
h = tfinal - t;
absh = abs(h);
done = true;
  end

  % LOOP FOR ADVANCING ONE STEP.
  nofailed = true;                   % no failed attempts
  while true
hA = h * A;
hB = h * B;
f(:,2) = feval(ode,t+hA(1),y+f*hB(:,1),args{:});
f(:,3) = feval(ode,t+hA(2),y+f*hB(:,2),args{:});
f(:,4) = feval(ode,t+hA(3),y+f*hB(:,3),args{:});
f(:,5) = feval(ode,t+hA(4),y+f*hB(:,4),args{:});
f(:,6) = feval(ode,t+hA(5),y+f*hB(:,5),args{:});

tnew = t + hA(6);
if done
   tnew = tfinal; % Hit end point exactly.
end
h = tnew - t;    % Purify h.

ynew = y + f*hB(:,6);
f(:,7) = feval(ode,tnew,ynew,args{:});
stats.nfevals = stats.nfevals + 6;

% Estimate the error.
if normcontrol
   normynew = norm(ynew);
   err = absh * (norm(f * E) / max(max(normy,normynew),threshold));
else
   err = absh * norm((f * E) ./ max(max(abs(y),abs(ynew)),threshold),inf);
end

% Accept the solution only if the weighted error is no more than the
% tolerance rtol.  Estimate an h that will yield an error of rtol on
% the next step or the next try at taking this step, as the case may be,
% and use 0.8 of this value to avoid failures.
if err > rtol                      % Failed step
   stats.nfailed = stats.nfailed + 1;
   if absh <= hmin
      msg = sprintf(['Failure at t=%e.  Unable to meet integration ' ...
                     'tolerances without reducing the step size below ' ...
                     'the smallest value allowed (%e) at time t.\n'],t,hmin);
      warning(msg);
      if haveoutfun
      feval(outfun,[],[],'done',outputArgs1{:});
      end
      if printstats                % print cost statistics
      fprintf('%g successful steps\n', stats.nsteps);
      fprintf('%g failed attempts\n', stats.nfailed);
      fprintf('%g function evaluations\n', stats.nfevals);
      fprintf('%g partial derivatives\n', stats.npds);
      fprintf('%g LU decompositions\n', stats.ndecomps);
      fprintf('%g solutions of linear systems\n', stats.nsolves);
      end
      if nargout > 0
      if soloutRequested
         sol.x = tout(1:nout).';
         sol.y = yout(1:nout,

.';
sol.solver = 'ode45';
sol.idata.f3d = f3d(1:nout,:,

;
         if haveeventfun
            sol.xe = teout.';
            sol.ye = yeout.';
            sol.ie = ieout.';
         end
         varargout{1} = sol;
      else
         varargout{1} = tout(1:nout);
         varargout{2} = yout(1:nout,

;
         if haveeventfun
            varargout{3} = teout;
            varargout{4} = yeout;
            varargout{5} = ieout;
            if ~FcnHandlesUsed
            varargout{6} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                              stats.npds; stats.ndecomps; stats.nsolves];
            end
         else
            if ~FcnHandlesUsed
            varargout{3} = [stats.nsteps; stats.nfailed; stats.nfevals; ...
                              stats.npds; stats.ndecomps; stats.nsolves];
            end
         end
      end
      end
      return;
   end

   if nofailed
      nofailed = false;
      absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
   else
      absh = max(hmin, 0.5 * absh);
   end
   h = tdir * absh;
   done = false;

else                               % Successful step
   break;

end
  end
  stats.nsteps = stats.nsteps + 1;

  if haveeventfun
[te,ye,ie,valt,stop] = ...
      odezero(@ntrp45,eventfun,eventargs,valt,t,y,tnew,ynew,t0,h,f);
nte = length(te);
if nte > 0
   if soloutRequested | nargout > 2
      teout = [teout; te];
      yeout = [yeout; ye.'];
      ieout = [ieout; ie];
   end
   if stop                         % stop on a terminal event
      if ~soloutRequested % preserve tnew,ynew when sol requested
      tnew = te(nte);
      ynew = ye(:,nte);
      end
      done = true;
   end
end
  end

  if nargout > 0
oldnout = nout;
if outflag == 3                   % computed points, with refinement
   nout = nout + refine;
   if nout > length(tout)
      tout = [tout; zeros(chunk,1)];  % requires chunk >= refine
      yout = [yout; zeros(chunk,neq)];
   end
   i = oldnout+1:nout-1;
   tout(i) = t + (tnew-t)*S;
   yout(i,

= ntrp45(tout(i),t,y,[],[],h,f).';
tout(nout) = tnew;
yout(nout,

= ynew.';
elseif outflag == 2                % computed points, no refinement
   nout = nout + 1;
   if nout > length(tout)
      tout = [tout; zeros(chunk,1)];
      yout = [yout; zeros(chunk,neq)];
      if soloutRequested
      f3d = [f3d; zeros(chunk,neq,7)];
      end
   end
   tout(nout) = tnew;
   yout(nout,

= ynew.';
if soloutRequested
f3d(nout,:,

= f;
   end
elseif outflag == 1                % output only at tspan points
   while next <= ntspan
      if tdir * (tnew - tspan(next)) < 0
      if haveeventfun & done
         nout = nout + 1;
         tout(nout) = tnew;
         yout(nout,

= ynew.';
      end
      break;
      elseif tnew == tspan(next)
      nout = nout + 1;
      tout(nout) = tnew;
      yout(nout,

= ynew.';
      next = next + 1;
      break;
      end
      nout = nout + 1;             % tout and yout are already allocated
      tout(nout) = tspan(next);
      yout(nout,

= ntrp45(tspan(next),t,y,[],[],h,f).';
      next = next + 1;
   end
end

if haveoutfun
   i = oldnout+1:nout;
   if ~isempty(i) & (feval(outfun,tout(i),yout(i,outputs).',outputArgs{:}) == 1)
      if soloutRequested
      sol.x = tout(1:nout).';
      sol.y = yout(1:nout,

.';
sol.solver = 'ode45';
sol.idata.f3d = f3d(1:nout,:,

;
      if haveeventfun
         sol.xe = teout.';
         sol.ye = yeout.';
         sol.ie = ieout.';
      end
      varargout{1} = sol;
      else
      varargout{1} = tout(1:nout);
      varargout{2} = yout(1:nout,

, ynew(outputs)],outputArgs{:}) == 1
      return;
   end
elseif outflag == 2
   if feval(outfun,tnew,ynew(outputs),outputArgs{:}) == 1
      return;
   end
elseif outflag == 1                % output only at tspan points
   ninterp = 0;
   while next <= ntspan
      if tdir * (tnew - tspan(next)) < 0
      if haveeventfun & done
         ninterp = ninterp + 1;
         tinterp(ninterp,1) = tnew;
         yinterp(:,ninterp) = ynew;
      end
      break;
      elseif tnew == tspan(next)
      ninterp = ninterp + 1;
      tinterp(ninterp,1) = tnew;
      yinterp(:,ninterp) = ynew;
      next = next + 1;
      break;
      end
      ninterp = ninterp + 1;
      tinterp(ninterp,1) = tspan(next);
      yinterp(:,ninterp) = ntrp45(tspan(next),t,y,[],[],h,f);
      next = next + 1;
   end
   if ninterp > 0
      if feval(outfun,tinterp(1:ninterp),yinterp(outputs,1:ninterp),outputArgs{:}) == 1
      return;
      end
   end
end
  end

  % If there were no failures compute a new h.
  if nofailed
% Note that absh may shrink by 0.8, and that err may be 0.
temp = 1.25*(err/rtol)^pow;
if temp > 0.2
   absh = absh / temp;
else
   absh = 5.0*absh;
end
  end

  % Advance the integration one step.
  t = tnew;
  y = ynew;
  if normcontrol
normy = normynew;
  end
  f(:,1) = f(:,7);          % Already evaluated feval(ode,tnew,ynew,args)

end

if haveoutfun
  feval(outfun,[],[],'done',outputArgs1{:});
end

if printstats                         % print cost statistics
  fprintf('%g successful steps\n', stats.nsteps);
  fprintf('%g failed attempts\n', stats.nfailed);
  fprintf('%g function evaluations\n', stats.nfevals);
  fprintf('%g partial derivatives\n', stats.npds);
  fprintf('%g LU decompositions\n', stats.ndecomps);
  fprintf('%g solutions of linear systems\n', stats.nsolves);
end

if nargout > 0
  if soloutRequested
sol.x = tout(1:nout).';
sol.y = yout(1:nout,

.';
sol.solver = 'ode45';
sol.idata.f3d = f3d(1:nout,:,

;
if haveeventfun
   sol.xe = teout.';
   sol.ye = yeout.';
   sol.ie = ieout.';
end
varargout{1} = sol;
  else
varargout{1} = tout(1:nout);
varargout{2} = yout(1:nout,