数学建模基础（Ⅱ）

我将数学建模基础分为三个部分浅谈自己的理解，同时分享MATLAB的演示和学习过程。第一部分为MATLAB基础内容，二三部分为基本数学方法的原理和实现，第二部分对数学原理有较为详细的探究。本部分包括以下三方面内容。

方程求根
数值积分
微分方程

---

方程求根

对于方程$f(x)=0$，当$f(x)$是线性函数时，方程就是线性方程，对于线性方程和线性方程组的求解，理论和数值求解方法丰富，在下文线性代数部分进行学习；当$f(x)$是非线性函数时，方程就是非线性方程，非线性方程的$f(x)$复杂多样，尚无一般的解析解法。比如$a_n{x}^n+a_{n-1}{x}^{n-1}+\cdot \cdot \cdot +a_{1}x+a_0=0$和$lnx-xsin\frac{\pi }{2}x+1=0$，前者为n次代数方程，$n\ge 5$时便没有求根公式可用，后者为超越方程，很难求得解析解。在无法求得解析解的时候，设计迭代方程，计算出方程满足精度要求的近似根，就可以认为在满足实际需求这一层面，计算问题已经解决。

一、迭代方法

迭代格式是，将方程$f(x)=0$转化为等价的方程$x=\phi (x)$，据此构造出$x_{n+1}=\phi (x_n)$，以此式迭代，选定一个任意的$x_0$，每个迭代式的$x_{n+1}$构成一个迭代数列 { x n } \{x_n\} {xn​}，若数列有极限$\underset{n\to \infty }{\lim }{x}_n=x^{\ast }$，则$x^{\ast }$即为方程的根。乍看可能有些抽象，试验一下便会理解。
简单迭代法
由迭代格式容易理解，简单迭代法即凑出符合迭代格式的简单式子，然后按照精度或者迭代次数进行迭代，不同的迭代式子会影响到结果的准确性和迭代次数，下文中提到的求根方法都是在寻找更准确高效的迭代式。
二分法
原理是，将区域二分，根据零点定理，判断根在哪个分段内，重复进行二分，直到满足精度ε得到近似解。算法的流程是，赋初值$a_k=a,b_k=b,k=0$即最开始的区间是$[a,b]$，令$x_k=\frac{a_k+b_k}{2}$，计算$f(x_k)$，然后根据零点定理确定二分后的区间，「若$f(x_k)=0$则$x_k$是$f(x)=0$的根，若$f(a_k)f(x_k)<0$则$a_{k+1}=a_k,b_{k+1}=x_k$，否则$a_{k+1}=x_k,b_{k+1}=b_k$。」令$k=k+1$，$x_k=\frac{a_k+b_k}{2}$，产生新的二分点，如果$|f(x_k)-0|\le \epsilon$退出计算，反之继续使用零点定理判断。根据原理和步骤，简单的计算一下$x^3+x-1=0$在$[0,1]$内的近似根，精度ε要求$10^{-5}$。
牛顿法
原理是，从迭代初值$x_0$对应找到方程$f(x)=0$中的函数$f(x)$上的一条切线，切线与$x$轴的交点成为下一个迭代$x$值，又对应函数上的一条切线，可以发现每次切线与$x$轴交点都会更逼近函数的零点，以此求得方程根。主要思想是，判断实根，「设函数$f(x)$在区间$[a,b]$上有二阶导数，$f(a)f(b)<0$，并且$f^{\prime }(x)、f^{\prime \prime }(x)$在区间内不变号，则方程$f(x)=0$在区间$(a,b)$内有且仅有一个实根$x^{\ast }$。」设迭代初值$x_0$，过$(x_0,f(x_0))$的切线方程为$y=f(x_0)+f^{\prime }(x_0)(x-x_0)$，令$y=0$可得切线与$x$轴交点为$x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$，则根据原理可以推出迭代式子即为$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$，这就是牛顿迭代公式。可证牛顿迭代式收敛。
弦截法
牛顿法每次需要计算$f(x)$和$f^{\prime }(x)$，遇到导数计算困难、在根附近导数值非常小、方程有多重根等问题时，有很大局限性，弦截法通过使用割线代替切线，因为$f^{\prime }(x_0)$零点可以直接得出，所以在$n\ge 1$时，让$f^{\prime }(x_n)=\frac{f(x_n)-f(x_0)}{x_n-x_0}$，在牛顿迭代式的基础上，可以得到单点弦截法迭代公式，$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}=x_n-\frac{f(x_n)(x_n-x_0)}{f(x_n)-f(x_0)}$，若将$x_0$改为$x_{n-1}$，即两个$x$值都更新，则会得到割线更新公式，$f^{\prime }(x_n)=\frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}$，两点弦截法迭代公式$x_{n+1}=x_n-\frac{f(x_n)(x_n-x_{n-1})}{f(x_n)-f(x_{n-1})}$。

二、算法加速

可以采用松弛法加速迭代收敛，或者将发散的迭代收敛。松弛法是指对迭代式$x_{n+1}=\phi (x_n)$选取松弛因子ω处理为$x_{n+1}=\omega \phi (x_n)+(1-\omega )x_n$。也可以通过Aitken法或Steffenson法加速，Aitken法避免了松弛法计算${\phi }^{\prime }(x_n)$的不便，利用迭代式$x_{n+1}=\phi (x_n)$方程根$x^{\ast }$，经过一系列计算得到复杂迭代式，Steffenson法也较复杂，不再详细计算。

三、MATLAB方程求解

MATLAB函数
$[x1,x2,\cdot \cdot \cdot ,xn]=solve{(}^{\prime }F{1}^{\prime }{,}^{\prime }F{2}^{\prime },\cdot \cdot \cdot {,}^{\prime }F{n}^{\prime }{,}^{\prime }var{1}^{\prime }{,}^{\prime }var{2}^{\prime },\cdot \cdot \cdot {,}^{\prime }var{n}^{\prime })$可以求解n个线性或非线性方程组成的方程组问题，F指各个方程，var指各个求解变量，x指对应var的求解结果。
另外$fzero$，$fsolve$，$roots$函数也都可用于求方程和方程组解。
图形法
绘制出函数曲线后，放大局部判断解的大体范围，以选择迭代初值。

四、代码

数学建模基础ⅠⅡⅢ中的绝大部分代码均来自《MATLAB与数学实验》一书，主要目的是练习，自己也会在书本代码之外写一些其他题目。

#Chapter2.m
syms x y y1 y2 y3 y4 a b c

%来一个线性方程，下面用各种迭代法求其解
fun=x^3+x-1;

vpa(exam2_1(fun,0,1,1e-5),7)

[y1,y2,y3,y4]=exam2_2(10)

vpa(exam2_3(fun,1,10),7)

vpa(exam2_4(fun,1e-5,100,1.8,1.9),7)

%直接用MATLAB函数求解
double(solve(fun,x))%solve解若含root需用double转换

f1=x^2+2*y-a;
f2=3*x-4*y-b;%解方程组a=x^2+2*y;b=3*x-4*y
[x,y]=solve(f1,f2,x,y)

%fzero和fsolve不能直接用上述符号方法建立的方程，可以用内联函数建立，也可以用@匿名函数，然后求非线性方程的最优解
f=inline('x+2*sin(x)*exp(x)-1');
fzero(f,[0,1]),fzero(@exam2_5,[0,1])

fsolve(@exam2_6,[0,0])

%作图以确定x^5-6x^4+5x^2-3x+8=0根的大致范围
x=-2:0.01:6;
y=x.^5-6*x.^4+5*x.^2-3*x+8;
figure,plot(x,y);
grid on;

figure,ezplot('x^5-6*x^4+5*x^2-3*x+8',[-2,6]);
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_1.m
function y=exam2_1(ff,m,n,er)
%迭代算法二分法算x^3+x-1=0在[0,1]的近似根（精度10^-5）
%fun=x^3+x-1，[m,n]内，精度为er。
syms x xk
a=m;b=n;k=0;
while b-a>er
    xk=(a+b)/2;
    fx=subs(ff,x,xk);
    fa=subs(ff,x,a);
    k=k+1;
    if fx==0
        y(k)=xk;
        break;
    elseif fa*fx<0
        b=xk;
    else
        a=xk;
    end
    y(k)=xk;
end
% figure,plot(y,'.-');
% grid on当然可以画迭代折线图
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_2.m
function [y1,y2,y3,y4]=exam2_2(n)
%简单迭代方式求x^2-2=0的根，n表示迭代次数。
y1=ones(1,n);y2=ones(1,n);y3=ones(1,n);y4=ones(1,n);
for i=1:n:1
    y1(i+1)=2/y1(i);
    y2(i+1)=0.5*(y2(i)+2/y2(i));
    y3(i+1)=(2*y3(i)-y3(i)^2+2)^(1/3);
    y4(i+1)=1+1/(1+y4(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_3.m
function y=exam2_3(ff,xx,n)
%牛顿迭代法计算x^3+x-1=0在[0,1]内的近似根
%无论哪种迭代方法的初值，都是越接近所求值越好。
syms x
fundiff=diff(ff);
xi=xx;
for i=1:1:10
    funA=subs(ff,x,xi);
    fundiffA=subs(fundiff,x,xi);
    xi=xi-funA/fundiffA;
end
y=xi;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_4.m
function y=exam2_4(ff,er,n,xa,xb)
%两点弦截法，输入分别为方程，精度，步数，迭代初值ab。
%y是向量，表示迭代过程中的各点，k时实际迭代步数。
syms x xk
x0=xa;x1=xb;
y(1)=xa;y(2)=xb;
k=2;
while abs(x1-x0)>er&&k<n
    fx0=subs(ff,x,x0);
    fx1=subs(ff,x,x1);
    x2=x1-fx1*(x1-x0)/(fx1-fx0);
    k=k+1;
    y(k)=x2;
    x0=x1;
    x1=x2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_5.m
function f=exam2_5(x)
f=x+2*sin(x)*exp(x)-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam2_6.m
function F=exam2_6(x)
F(1)=x(1)-sin(x(1))-5*x(2);
F(2)=2*x(1)+x(2)-cos(x(2));
end

数值积分

数值积分，也就是求定积分的近似值。定积分的定义是高数的基本内容，用微分和极限的思想将曲线和坐标系包围成的区域趋近无穷地划分为一个个无限微小的矩形，分别计算后求和得到定积分。如果被积函数f(x)有初等函数形式的原函数F(x)，那么对于定积分${\int }_a^{b}f(x)dx$使用Newto-Leibniz Formula可以得到公式${\int }_a^{b}f(x)dx=F(x){|}_a^b=F(b)-F(a)$，所谓的不定积分也就是原函数，作为一个函数集存在的F(x)，与定积分的关系从此公式也可以得出；但若被积函数f(x)不具有初等函数形式的原函数，或其理论存在但不能具体确定时，则考虑数值计算定积分的方法；另外，有时被积函数f(x)都不会明确给出解析表达式，而是以表格数据、曲线图形的形式给出，这时也只能考虑数值积分的方法计算相应定积分。

一、数值积分方法

一重积分
从上文积分的定义可以得到一个微分求和的式子${\int }_a^{b}f(x)dx=\underset{n\to \infty }{\lim }\frac{1}{n}(b-a)\sum _{i=1}^{n}f({\xi }_i)$，其中$\frac{1}{n}(b-a)$是每个近似小矩形的宽，${\xi }_i$由x决定，决定了矩形（其他近似图形也与之类似）的长，最终的近似结果也由${\xi }_i$起决定性作用，按照${\xi }_i$的不同进行方法分类，分别是左矩形法${\xi }_i=x_i$、右矩形法${\xi }_i=x_{i-1}$、中矩形法${\xi }_i=\frac{x_{i-1}+x_i}{2}$、梯形法$f({\xi }_i)=\frac{f(x_{i-1})+f(x_i)}{2}$、辛普森法$f({\xi }_i)=\frac{f(x_{i-1})+4f(\frac{x_{i-1}+x_i}{2})+f(x_i)}{6}$。
二重积分
二重积分从一元函数变为二元函数，也就是在三维坐标系中求形如z=f(x,y)函数的积分，需要对x和y两个变量求积，从原来求一个由曲线和坐标轴构成的曲边图形面积变成了求以函数面为顶的曲顶柱体的体积，需要求积分的微分单元就从矩形变成了长方体（或其他近似体），将XOY面上z(x,y)的投影记作D，将D划分为无限个小区域记为σ，每一个小柱体的体积就为$f(x_i,y_i)\Delta {\sigma }_i$，那么二重积分公式即为${\iint }_{D}f(x,y)d\sigma =\underset{\lambda \to 0}{\lim }\sum _{i=1}^{n}f(x_i,y_i)\Delta {\sigma }_i$，λ即每一块σ的最大值。求二重积分可以使用单点柱体法，就把小区域σ近似为等分的无穷多个矩形，并取$f(x_i,y_i)\Delta {\sigma }_i$中的$(x_i,y_i)$为${\sigma }_i$的中心点，公式就变成了${\iint }_{D}f(x,y)d\sigma =\underset{\lambda \to 0}{\lim }\sum _{i=1}^{n}f(x_i,y_i)\Delta x_i\Delta y_i$。

二、MATLAB数值积分

通过MATLAB内置函数求积分。
一重积分
$trapz(X,Y)$表示使用梯形方法计算Y对X的积分值，X和Y是维度相同的向量。
$quad(fun,a,b,tol,trace)，quadl$两个函数是相同的形式，分别采用辛普森方法、递推自适应Lobatto法求数值积分，fun为被积函数，a、b为积分上下限，tol为精度，默认时为1.0e-6。
$[q,errbnd]=quadgk(fun,a,b,parami,vari,...)$此函数采用自适应Gauss-Kronrod方法进行数值积分计算，适用于高精度积分、震荡数值积分，支持无穷区间，并且可以处理端点包含奇点的情况，同时还支持沿着不连续函数积分，复数域线性路径的围道积分法。其中fun表示被积函数a、b表示积分限，param、val为函数的其他控制参数，包括绝对误差、相对误差、路径点、最大积分区间数等。此外，q、errbnd分别是积分值和误差。
二重积分
$int的嵌套$int函数可以求一重积分，嵌套多层便可以简单求多重积分。
$dblquad(fun,xmin,xmax,ymin,ymax,tol,method)$表示fun在xmin<=x<=xmax，ymin<=y<=ymax范围内求解定积分，tol为误差，默认值也为1.0e-6，method为数值积分方法，例如quadl，默认为quad。

三、 代码

#Chapter3.m
syms x y c

%一重的
fun1=1/(1+x^2);
[ss,er]=exam3_1(fun1,1,0,1,50);
vpa([ss,er],7)

%二重的
[ss,er]=exam3_2(0,1,0,1,100,100);
vpa([ss,er],7)

[ss,er]=exam3_3(0,1,500,0.1);
vpa([ss,er],7)

%直接用MATLAB函数求解
%一重的
dx=pi/1000;
xx=0:dx:2*pi;
fun2=exam3_4(x);%exp(-0.5*x).*sin(x+pi/4)
y=subs(fun2,x,xx);%不用subs直接代入可以少写一步
trapz(y)*dx

fun3=exam3_5(x);%4./(x.^2+1)
quad(@exam3_5,0,1)

fun4=exam3_6(x,c);%1./(x.^3-2*x-c)
quadl(@exam3_6,0,2,[],[],10)

fun5=exam3_7(x);%exp(-x.^2).*(log(x)).^2
quadgk(@exam3_7,0,inf)

%二重的
syms x y
fun6=x*y+exp(x);%x*y+exp(x)
int(int(fun6,x,0,1),y,0,2)

fun7=exam3_8(x,y);%x.*exp(y)+y.*sin(x)
dblquad(@exam3_8,0,pi,0,1,1e-8,@quadl)%默认@quad
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_1.m
function [ss,er]=exam3_1(ff,t,a,b,n)
%t控制积分方法；[a,b]为积分区间；n是分割的份数，影响的是准确度。
%ss积分近似值，er误差。
dx=(b-a)/n;
syms x
xx=a:dx:b;
ss=0;
if t==1%左矩阵法
    for i=1:n
        f1=subs(ff,x,xx(i));
        ss=ss+dx*f1;
    end
elseif t==2%右矩阵法
    for i=1:n
        f1=subs(ff,x,xx(i+1));
        ss=ss+dx*f1;
    end
elseif t==3%中矩阵法
    for i=1:n
        f1=subs(ff,x,(xx(i)+xx(i+1))/2);
        ss=ss+dx*f1;
    end
elseif t==4%梯形法
    for i=1:n
        f1=subs(ff,x,xx(i));
        f2=subs(ff,x,xx(i+1));
        ss=ss+dx*(f1+f2)/2;
    end
else%辛普森法
    for i=1:n
        f1=subs(ff,x,xx(i));
        f2=subs(ff,x,xx(i+1));
        f3=subs(ff,x,(xx(i)+xx(i+1))/2);
        ss=ss+dx*(f1+f2+4*f3)/6;
    end
end
er=abs(ss-pi/4);%误差就是算出1/(1+x^2)积分理论解与之相比
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_2.m
function [ss,er]=exam3_2(a,b,c,d,m,n)
%矩形域单点柱体法计算定积分x*y+y^3，1<=x<=1，0<=y<=1
%[a,b]为积分区域x范围，[c,d]为积分区域y范围，m、n表示分割份数
%ss积分近似值，er误差
dx=(b-a)/m;
dy=(d-c)/n;
x=a:dx:b;
y=c:dy:d;
ss=0;
for i=1:1:m
    xx=(x(i)+x(i+1))/2;
    for j=1:1:n
        yy=(y(j)+y(j+1))/2;
        fxy=xx*yy+yy^3;%被积函数，这里换成subs，将符号函数作为输入参数计算就速度极慢，不知道该怎么写。
        ss=ss+fxy*dx*dy;
    end
end
er=abs(ss-0.5);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_3.m
function [ss,er]=exam3_3(a,b,n,k)
%X型或Y型二重积分计算xy在y=x与y=x^2围成区域内的积分
%[a,b]为定值范围，n表示分割份数，k表示dy=k*dx
%ss积分近似值，er误差
dx=(b-a)/n;
dy=k*dx;
x=a:dx:b;
ss=0;
for i=1:1:n
    xx=(x(i)+x(i+1))/2;
    y1=x(i)^2;
    y2=x(i);
    yk=y1:dy:y2;
    p=length(yk);
    if p>1
        for j=1:1:p-1
            yy=(yk(j)+yk(j+1))/2;
            fxy=xx*yy;%被积函数
            ss=ss+fxy*dx*dy;
        end
    end
end
er=abs(ss-1/24)*24;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_4.m
function y=exam3_4(x)
y=exp(-0.5*x).*sin(x+pi/4);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_5.m
function y=exam3_5(x)
y=4./(x.^2+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_6.m
function y=exam3_6(x,c)
y=1./(x.^3-2*x-c);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_7.m
function y=exam3_7(x)
y=exp(-x.^2).*(log(x)).^2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam3_8.m
function f=exam3_8(x,y)
f=x.*exp(y)+y.*sin(x);
end

微分方程

本部分内容学习自《数值分析》和《高等数学讲义》，大学开设的「数值计算」课程实在是水，只掌握了数值计算和分析方法的大概轮廓，一学期下来除了掌握了MATLAB之外几乎一无所获，如今自行学习，希望能学有所得。
包括建模题的很多实际问题中，需要考虑变化率问题，往往引入导数（或微分）来建立相关的数学模型，进行理论分析和数值求解。含有未知函数的导数（或微分）的方程叫做微分方程，其中包含未知函数的（导数或微分的）最高阶数称为微分方程的阶，如果微分方程所有导数只对一个变量求导，就叫做常微分方程，否则为偏微分方程。n阶常微分方程一般形式为$F(x,y,y^{\prime },y^{\prime \prime },...,y^{(n)})=0$，如果y是n阶可导函数$y=y(x)$且满足上式$F(x,y(x),y^{\prime }(x),y^{\prime \prime }(x),...,y^{(n)}(x))=0$「有且只有这一个函数满足上式」则，称函数y=y(x)为微分方程的解。若微分方程的解中所含的常数个数等于方程阶数，则称其为微分方程的通解；当通解的常数C「所谓常数C指的是还可取任意值但还未确定的常数」都取了定值，得到的就是微分方程的特解。微分方程的初始条件就是用来确定通解里面常数的东西，比如落体运动的物理模型公式，$v{t}^2-v_0{t}^2=2gh$，积分过程中出现的$C_0$和$C_1$也就是初速度$v_0$和起始位置$h$即为初始条件，带有初始条件求$y=f(x)$所有初值的问题称为初值问题。n阶微分方程的初值问题需要n个初始条件，比如二阶常微分方程初值问题的格式为$\{\begin{array}{l}{y}^{\prime \prime }=f(x,y)\\ y^{\prime }(x_0)=y_0\\ y(x_0)=y_1\end{array}$。由于解析解难以求得或不存在，所以一般通过下述方法来解常微分方程初值问题的数值解，在精度允许范围内求出近似解。
另外求常微分方程数值解是指在求解区间内一系列离散点处给出真解的近似值，并不是拟合出方程的特解。

一、求常微分方程数值解的方法

设一阶常微分方程的初值问题为$\{\begin{array}{l}{y}^{\prime }=f(t,y),a\le t\le b\\ y(a)=\eta \end{array}$，首先要求函数f(x,y)满足Lipschitz条件，然后我们总会发现，在实际问题中所归结得到的初值问题往往不能用解析方法求解，求初值问题的解的一类数值方法是离散变量法，其过程就叫做离散化过程，而解初值问题得离散化方法通常有三种：差商代替导数、Taylor级数、数值积分，这里只讲差商代替导数方法，包括欧拉方法和Runge-Kutta法。
离散变量法其目的是将一个连续型问题化为一个离散型问题，实际上是一种思想而不是方法，比如求初值问题，先求初值问题的精确解$y=y(t)$在一系列的离散点$t_1,t_2,...,t_N$上的近似值有$y_1,y_2,...,y_N$，用$y(t_n)$表示$y(t)$在$t=t_n$处的值，则取$y(t_n)\simeq y_n,n=1,2,...,N$相邻两个节点间距$h_n=t_{n+1}-t_n,「h_n$通常$=h,t_0=a」$h称为步长为常数，也就是说离散点$t_0,t_1,t_2,...,t_N$是等距的，这个时候就有了$t_n=a+nh,n=0,1,...,N,Nh=b-a$，其实理解起来很简单，整个过程可以理解为把自变量等分为一个个微小的点，对应每个点的因变量取近似值，这样求解。
差商有递归定义其中关于节点的0阶差商定义为其函数值，，可以这样理解，微商「即导数」，是微分之商，微分即$d$，是变量的无穷小的改变量，导数就是两个变量的无穷小的改变量的比值。同理，差商，是差分之商，差分即$\Delta$，是变量的改变量，差商就是两个变量的改变量的比值。粗略的理解，两者的差别在于改变量的尺度不一样。改变量$\Delta$趋于零的时候，差商就成了导数。（语出自知乎李晓）一个k阶常微分方程转化为k阶差分方程的推导：上文初值问题$\{\begin{array}{l}{y}^{\prime }=f(t,y),a\le t\le b\\ y(a)=\eta \end{array}$中的初始条件$y(a)=y(t_0)=\eta$是给定的，因此可以计算出$y^{\prime }(t_0)=f(t_0,y(t_0))$，之后使用差商代替导数的步骤如下。

欧拉方法
欧拉法的原理是在小区间$[x_n,x_{n+1}]$上使用差商代替导数$y^{\prime }$，即$y^{\prime }\approx \frac{y(x_{n+1})-y(x_n)}{x_{n+1}-x_n}$，也就是上文说的用差分代替微分，用差商代替导数，对离散点$x_n$采用等距节点$h=x_n-x_{n-1},h为常数(x_n=x_0+nh,n=1,2,...以下同)$并假设$f(x,y)\approx f(x_n,y(x_n))$，由上述各式可以推导出$y(x_{n+1})\approx y(x_n)+hf(x_n,y(x_n))$，假定$y(x_n)\approx y_n$则可以得到$y(x_{n+1})$的近似计算方法：$y_{n+1}=y(x_n)+hf(x_n,y(x_n))$，其他推导过程见下图，引自《MATLAB与数学实验》。欧拉法又可以叫做欧拉折线法，是通过折线逼近曲线，也就是差商代替导数的一种方式，将问题离散化，并改变离散点的距离尺度，近似求解。

Runge-Kutta法
这里就不一点点敲LaTeX了，下图引自《MATLAB与数学实验》。

二、MATLAB求解微分方程

解析解
$dsolve{(}^{\prime }eqn{1}^{\prime }{,}^{\prime }eqn{2}^{\prime },...{,}^{\prime }con{1}^{\prime }{,}^{\prime }con{2}^{\prime },...{,}^{\prime }{v}^{\prime })$可以用来求解微分方程解析解，eqn1、eqn2、…是输入的方程组；con1、con2、…是初始条件；v表示关于其求导的变量，因为是常微分方程所以只有一个变量，省略时系统默认变量为t，在输入方程和初始条件时，要用Dny表示y关于自变量的n阶导数。
数值解
$[T,Y]=solver(odefun,tspan,y_0,options)$可以用来求解微分方程的数值解，options用于对求解算法和控制条件进行进一步设置。odefun为微分方程中的f(t,y)项；tspan为求解区间，要获得问题在指定点的解可以取$pan=[t_0,t_1,...,t_f]$，单调；y0为初始条件。solver为数值求解微分方程函数，也就是下列ode系列函数。

三、代码

#Chapter4.m
syms x y1 y2 y3 y4

%欧拉法求解常微分方程
[y1,y2,y3,y4]=exam4_1(0,1,100)

%Runge-Kutta法求解常微分方程
[y1,y2,y3]=exam4_2(0,1,100)

%直接用MATLAB函数求解
%解析解
%微分方程通解
dsolve('D1y=exp(2*t)-sin(3*t)')
dsolve('D1y=x^2+y','x')
%微分方程特解
dsolve('D1y=-y+2*x+1','y(0)=1','x')
%高阶微分方程
dsolve('D2y-5*Dy+6*y=x*exp(2*x)','x')
%微分方程组通解
[Y,Z]=dsolve('Dy=x+2*y-3*z','Dz=3*x+y-2*z','x')
%微分方程组特解
[x,y]=dsolve('D2x+Dy+3*x=cos(2*t)','D2y-4*Dx+3*y=sin(2*t)','Dx(0)=1/5','x(0)=0','Dy(0)=6/5','y(0)=0','t')

%数值解
[x,y]=ode45(@exam4_3,[0,0.5],1)
plot(x,y,'o-')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam4_1.m
function [y1,y2,y3,y4]=exam4_1(a,b,k)
%a是起始点，b是迭代终止点，h=(b-a)/k
h=(b-a)/k;x=a:h:b;[m,n]=size(x);
y1=zeros(m,n),y2=zeros(m,n),y3=zeros(m,n);
y1(1)=1;y2(1)=1; y3(1)=1;
for i=1:1:n-1%向前欧拉法
    x1=x(i);x2=x(i+1);
    y1(i+1)=y1(i)+h*(2*y1(i)/(x1+1)+(x1+1)^(5/2));%直接用f(x,y)，其他方法需要推导
end
for i=1:1:n-1%向后欧拉法
    x1=x(i);x2=x(i+1);
    yy=y2(i)+h*(2*y2(i)/(x1+1)+(x1+1)^(5/2));
    y2(i+1)=y2(i)+h*(2*yy/(x2+1)+(x2+1)^(5/2));
end
for i=1:1:n-1%改进欧拉法
    x1=x(i);x2=x(i+1);
    z1=2*y3(i)/(x1+1)+(x1+1)^(5/2);
    yy=y3(i)+h*z1;
    z2=2*yy/(x2+1)+(x2+1)^(5/2);
    y3(i+1)=y3(i)+h*0.5*(z1+z2);
end
for i=1:1:n%精确解
    y4(i)=(x(i)+1)^2*(2/3*(x(i)+1)^(3/2)+1/3);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam4_2.m
function [y1,y2,y3]=exam4_2(a,b,k)
%a是起始点，b是迭代终止点，h=(b-a)/k
h=(b-a)/k;x=a:h:b;[m,n]=size(x);
y1=zeros(m,n),y2=zeros(m,n),y3=zeros(m,n);
y1(1)=1;y2(1)=1; y3(1)=1;
for i=1:1:n-1%改进欧拉法
    x1=x(i);x2=x(i+1);
    z1=y1(i)-2*x1/y1(i);
    yy=y1(i)+h*z1;
    z2=yy-2*x2/yy;
    y1(i+1)=y1(i)+h*0.5*(z1+z2);
end
for i=1:1:n-1%四阶Runge-Kutta法
    x1=x(i);x2=x(i+1);xx=x1+h/2;
    k1=y2(i)-2*x1/y2(i);
    yy=y2(i)+h*k1/2;
    k2=yy-2*xx/yy;
    yy=y2(i)+h*k2/2;
    k3=yy-2*xx/yy;
    yy=y2(i)+h*k3;
    k4=yy-2*x2/yy;
    y2(i+1)=y2(i)+h/6*(k1+2*k2+2*k3+k4);
end
for i=1:1:n%精确解
    y3(i)=sqrt(1+2*x(i));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%exam4_3.m
function f=exam4_3(x,y)
f=-2*y+2*x^2+2*x;
end

总结

最近半年有些迷失了自我，推免这条路不知道还能不能走通，凭现在的GPA想有好出路简直是痴人说梦，拿国家级奖项、发表论文、完成省级立项全部顺利的情况下或许还有一丝希望，必须靠最后这一年刷绩点，考研这条路又有多艰难？大一的恶果终于自食，看不到未来。