## MIT-微分方程-matlab

### 欧拉方法ODE1

(a)   例子
F=@(t,y).06*y-t        %这里的.06就是0.06，我们输入t和y值，然后按公式计算

(b)   例子(logisctic equation)
a=2;
b=3;
F=@(t, y)a*y-b*y^2

function yout = ode1(F,t0,h,tfinal,y0)
% ODE1  A simple ODE solver.
%   yout = ODE1(F,t0,h,tfinal,y0) uses Euler's
%   method with fixed step size h on the interval
%      t0 <= t <= tfinal
%   to solve
%      dy/dt = F(t,y)
%   with y(t0) = y0.

%   Copyright 2014 - 2015 The MathWorks, Inc.

y = y0;
yout = y;
for t = t0 : h : tfinal-h
s = F(t,y);
y = y + h*s;
yout = [yout; y];  %每次得到的新的y都放在yout后面拼接起来
end

F=@(t,y) 2*y;
t0=0;
h=1;
tfinal=3;
y0=10;
ode1(F,t0,h,tfinal,y0)

F=@(t,y) .06*y;  %
t0=0;
h=1/12;
tfinal=10;
y0=100;
compound=ode1(F,t0,h,tfinal,y0)
format bank  %取两位小数，分别对应美元和美分
t=(0:h:tfinal)';   %转置一下
simple=y0+100*0.06*t;
plot(t,[compound,simple],'.')  %对比两种存钱方法

### 中间点法ODE2

ODE1的方法是用个所在点的斜率，推导下一个点的高度位置；ODE2是用相邻两个点中间的点的斜率确定下一个函数点的位置。方法为\begin{aligned} &s_{1}=f\left(t_{n}, y_{n}\right)\\ &s_{2}=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} s_{1}\right)\\ &y_{n+1}=y_{n}+h s_{2} \end{aligned}

function yout = ode2(F,t0,h,tfinal,y0)
% ODE2  A simple ODE solver.
%   yout = ODE2(F,t0,h,tfinal,y0) uses a midpoint
%   rule with fixed step size h on the interval
%      t0 <= t <= tfinal
%   to solve
%      dy/dt = F(t,y)
%   with y(t0) = y0.

%   Copyright 2014 - 2015 The MathWorks, Inc.

y = y0;
yout = y;
for t = t0 : h : tfinal-h
s1 = F(t,y);
s2 = F(t+h/2, y+h*s1/2);
y = y + h*s2;
yout = [yout; y];
end

F=@(t,y) sqrt(1-y^2);
t0=0;
h=pi/32;
tfinal=pi/2;
y0=0;
y=ode2(F,t0,h,tfinal,y0);
t=(0:h:tfinal)';   %转置一下
simple=sin(t);
plot(t,[y,simple],'o-')
axis([0 pi/2 0 1])
set(gca,'xtick',[0:pi/8:pi/2],'xticklabel',{'0','\pi/8','\pi/4','3\pi/8','\pi/2'})
xlabel('t'),ylabel('y')

### 经典龙格-库塔法ODE4

function yout = ode4(F,t0,h,tfinal,y0)
% ODE4  Classical Runge-Kutta ODE solver.
%   yout = ODE4(F,t0,h,tfinal,y0) uses the classical
%   Runge-Kutta method with fixed step size h on the interval
%      t0 <= t <= tfinal
%   to solve
%      dy/dt = F(t,y)
%   with y(t0) = y0.

%   Copyright 2014 - 2015 The MathWorks, Inc.

y = y0;
yout = y;
for t = t0 : h : tfinal-h
s1 = F(t,y);
s2 = F(t+h/2, y+h*s1/2);
s3 = F(t+h/2, y+h*s2/2);
s4 = F(t+h, y+h*s3);
y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
yout = [yout; y];
end

F=@(t,y) y^2-y^3;
t0=0;
y0=0.01;
tfinal=2/y0;
h=0.4;
y=ode4(F,t0,h,tfinal,y0);
t=(0:h:tfinal)';   %转置一下
plot(t,y,'.')

If the differential equation$$\frac{d y}{d t}=f(t)$$only involves $$t,$$ not $$y,$$ then its solution is just an integral,$$y(t)=\int f(t) d t$$In this situation, the Runge-Kutta method becomes a classic method of numerical quadrature. It should be familiar if you have ever studied such methods. Do you recognize it?

What is the solution to$$\frac{d y}{d t}=1+y^{2}$$with$$y(0)=0$$What happens if you try to use ode4 to solve this equation on the interval$$0 \leq t \leq 2$$

What happens with ode4 if the length of the interval, tfinal-to, is not evenly divisible by the step size, h? For example, if t0=0, tfinal =p i, and h=0 .1. This is one of the hazards of a fixed step size.

$$o(x^p)$$表示的是$$x^p$$的高阶无穷小，也就是当$$p$$趋近于零时，$$o(x^p)$$比$$x^p$$更快趋近于零。比如$$x^2$$是$$x$$的高阶无穷小，可以表示为$$x^2=o(x)$$。回忆一下什么叫作“同阶无穷小”，什么叫作“等价无穷小”(比如$$\sin x$$和$$x$$)。

对于经典的四阶龙格-库塔法，每一步的local error(局部截断误差，或者说每一步的误差)为$$h^5$$阶，总的global error为$$h^4$$阶，如果让步长减半为$$\frac { h}{ 2}$$，那么误差将会减少$$2^4$$。对于ODE1欧拉方法是一个一阶方法，意味着其局部截断误差（每步误差）正比于步长的平方，并且其全局截断误差正比于步长。

Modify order x to do further experiments involving the order of our ODE solvers.

### ODE23(估计误差，自动调整步长)

modern numercial method 不同于前面的ODE1欧拉方法或者四阶龙格-库塔方法，它的特点是：
(1)    Error estimate and step size control;
(2)   Fully accurate continuous interpolant(插值)

Bogacki-Shampine Order 3(2)    (ODE23的basis)\begin{aligned} s_{1} &=f\left(t_{n}, y_{n}\right) \\ s_{2} &=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} s_{1}\right) \\ s_{3} &=f\left(t_{n}+\frac{3}{4} h, y_{n}+\frac{3}{4} h s_{2}\right) \\ t_{n+1} &=t_{n}+h \\ y_{n+1} &=y_{n}+\frac{h}{9}\left(2 s_{1}+3 s_{2}+4 s_{3}\right) \\ s_{4} &=f\left(t_{n+1}, y_{n+1}\right) \\ e_{n+1} &=\frac{h}{72}\left(-5 s_{1}+6 s_{2}+8 s_{3}-9 s_{4}\right) \end{aligned}每一步计算的误差都去和这里的$$e_{n+1}$$进行比较，如果实际误差比它大，那么减小步长(多插入点)。

function [tout,yout] = ode23tx(F,tspan,y0,arg4,varargin)
%ODE23TX  Solve non-stiff differential equations.  Textbook version of ODE23.
%
%   ODE23TX(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system
%   of differential equations dy/dt = f(t,y) from t = T0 to t = TFINAL.
%   The initial condition is y(T0) = Y0.
%
%   The first argument, F, is a function handle or an anonymous function
%   that defines f(t,y).  This function must have two input arguments,
%   t and y, and must return a column vector of the derivatives, dy/dt.
%
%   With two output arguments, [T,Y] = ODE23TX(...) returns a column
%   vector T and an array Y where Y(:,k) is the solution at T(k).
%
%   With no output arguments, ODE23TX plots the emerging solution.
%
%   ODE23TX(F,TSPAN,Y0,RTOL) uses the relative error tolerance RTOL
%   instead of the default 1.e-3.
%
%   ODE23TX(F,TSPAN,Y0,OPTS) where OPTS = ODESET('reltol',RTOL, ...
%   'abstol',ATOL,'outputfcn',@PLOTFUN) uses relative error RTOL instead
%   of 1.e-3, absolute error ATOL instead of 1.e-6, and calls PLOTFUN
%   instead of ODEPLOT after each successful step.
%
%   More than four input arguments, ODE23TX(F,TSPAN,Y0,RTOL,P1,P2,...),
%   are passed on to F, F(T,Y,P1,P2,...).
%
%   ODE23TX uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23).
%
%   Example
%      tspan = [0 2*pi];
%      y0 = [1 0]';
%      F = @(t,y) [0 1; -1 0]*y;
%      ode23tx(F,tspan,y0);
%

%   Copyright 2012 - 2015 The MathWorks, Inc.

% Initialize variables.

rtol = 1.e-3;
atol = 1.e-6;
plotfun = @odeplot;
if nargin >= 4 & isnumeric(arg4)
rtol = arg4;
elseif nargin >= 4 & isstruct(arg4)
if ~isempty(arg4.RelTol), rtol = arg4.RelTol; end
if ~isempty(arg4.AbsTol), atol = arg4.AbsTol; end
if ~isempty(arg4.OutputFcn), plotfun = arg4.OutputFcn; end
end
t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
plotit = (nargout == 0);
threshold = atol / rtol;
hmax = abs(0.1*(tfinal-t0));
t = t0;
y = y0(:);

% Initialize output.

if plotit
plotfun(tspan,y,'init');
else
tout = t;
yout = y.';
end

% Compute initial step size.

s1 = F(t, y, varargin{:});
r = norm(s1./max(abs(y),threshold),inf) + realmin;
h = tdir*0.8*rtol^(1/3)/r;

% The main loop.

while t ~= tfinal

hmin = 16*eps*abs(t);
if abs(h) > hmax, h = tdir*hmax; end
if abs(h) < hmin, h = tdir*hmin; end

% Stretch the step if t is close to tfinal.

if 1.1*abs(h) >= abs(tfinal - t)
h = tfinal - t;
end

% Attempt a step.

s2 = F(t+h/2, y+h/2*s1, varargin{:});
s3 = F(t+3*h/4, y+3*h/4*s2, varargin{:});
tnew = t + h;
ynew = y + h*(2*s1 + 3*s2 + 4*s3)/9;
s4 = F(tnew, ynew, varargin{:});

% Estimate the error.

e = h*(-5*s1 + 6*s2 + 8*s3 - 9*s4)/72;
err = norm(e./max(max(abs(y),abs(ynew)),threshold),inf) + realmin;

% Accept the solution if the estimated error is less than the tolerance.

if err <= rtol
t = tnew;
y = ynew;
if plotit
if plotfun(t,y,'');
break
end
else
tout(end+1,1) = t;
yout(end+1,:) = y.';
end
s1 = s4;     % Reuse final function value to start new step.
end

% Compute a new step size.

h = h*min(5,0.8*(rtol/err)^(1/3));

% Exit early if step size is too small.

if abs(h) <= hmin
warning('Step size %e too small at t = %e.\n',h,t);
t = tfinal;
end
end

if plotit
plotfun([],[],'done');
end

F=@(t,y)y;      %微分方程，导函数等于函数值本身(其实是指数函数)
ode23(F,[0 1],1)      %  interval from 0 to 1, 初始值为1，没有output argument(输出参数)

F=@(t,y)y;
[t,y]=ode23(F,[0 1],1)   %对比上一个例子加入输出参数
plot(t,y,'o-')                      %的到的图和上面一样

a=0.25;
y0=15.9;
F=@(t,y)2*(a-t)*y^2
[t,y]=ode23(F,[0 1],y0)
plot(t,y,'-',t,y,'.','markersize',24)

h=diff(t);
plot(t(2:end),h,'.','markersize',24)这两行代码，可以得到如下图像