文章目录

一、数值微分与数值积分
- 1、数值微分
- 2、数值积分
二、线性方程求解
- 1、直接法
- 2、迭代解法
- 3、直接法与迭代法的对比
三、非线性方程求解与函数极值计算
- 1、非线性方程数值求解
- 2、函数极值的计算
四、常微分方程数值求解
- 1、常微分方程数值求解的一般概念
- 2、常微分方程数值求解函数
- 3、刚性问题
五、常微分方程应用举例
- 1、Lotka-Volterra模型
- 2、Lotka-Volterra改进模型
总结

一、数值微分与数值积分

1、数值微分

MATLAB提供了求向前差分的函数diff，调用格式：
(1)dx=diff(x)：计算向量x的向前差分，dx(i)=x(i+1)-x(i),i=1,2,……,n-1。
(2)dx=diff(x,n)：计算向量x的n阶向前差分，例如，diff(x,2)=diff(diff(x))。
(3)dx=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(默认状态)，按列计算差分；dim=2时，按行计算差分。
注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样对于矩阵，差分后的矩阵比原矩阵少了一行或者一列。
另外，计算差分之后，可以用f(x)在某点的差商作为其导数的近似值。
例子：设f(x)=sinx，在[0,2Π]范围内随机采样，计算f’(x)的近似值，并与理论f’(x)=cosx进行比较

x=linspace(0,2*pi,5000);
y=sin(x);
f1=diff(y)./diff(x);
f2=cos(x(1:end-1));%差分向量元素比原向量元素少了一个
plot(x(1:end-1),f1,'r',x(1:end-1),f2,'b');
d=norm(f1-f2)%求f1,f2的2范数

2、数值积分

(1)数学原理
常用的求定积分的方式是利用牛顿-莱布尼兹公式：
∫abf(x)dx\int_{a}^{b}f(x)dx∫abf(x)dx=F(b)-F(a)，其中F(x)为f(x)的原函数。
但是当被积函数的原函数无法用初等函数表示，或者被积函数是用离散的形式给出，这个时候需要用数值解法求解定积分。
求定积分的数值解法很多种，如梯形法、辛普森法、高斯求积法等等。基本思想都是将积分区间分成n个子区间，这样定积分问题就分解成每个子区间上积分后求和的问题：
S=∫abf(x)dx\int_{a}^{b}f(x)dx∫abf(x)dx=∑i=1n\sum_{i=1}^n∑i=1n∫xixi−1f(x)dx\int_{x_i}^{x_{i-1}}f(x)dx∫xixi−1f(x)dx
(2)数值积分的实现

基于自适应辛普森方法
[I,n]=quad(filename,a,b,tol,trace)
基于自适应Gauss-Lobattofangf
[I,n]=quadl(filename,a,b,tol,trace)
其中，filename是被积函数名；a和b分别是定积分的下限和上限，积分限[a,b]必须是有限的，不能为无穷大；tol用来控制积分精度，默认取tol=10−610^{-6}10−6；trace控制是否展现积分过程，若取非0则展现，取0则不展现，默认取0；返回参数I即定积分的值，n为被积函数的调用次数。
例子：分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下比较被积函数的调用次数。

format long
f=@(x)4./(1+x.^2);
[I1,n1]=quad(f,0,1,1e-8)
[I2,n2]=quadl(f,0,1,1e-8)
(atan(1)-atan(0))*4%输出理论值

基于全局自适应积分方法
I=integral(filename,a,b)
其中，I是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。
例子：求定积分I=∫1∞1x1−ln⁡2xdx\int_{1}^{∞}{\frac 1{x\sqrt{1-\ln^2x}}}dx∫1∞x1−ln2x1dx
基于自适应的Gauss-Lobattofangf(高斯-克朗罗德)方法
[I,err]=quadgk(filename,a,b)
其中，err为返回近似误差范围，其他参数与quad函数相同。积分上下限可以是无穷大(-Inf或Inf)，也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。
例子：求定积分∫2π+∞1x2sin⁡1xdx\int_{\frac 2\pi}^{+∞}{{1 \over x^2}{\sin{1\over x}}}dx∫π2+∞x21sinx1dx

f=@(x)sin(1./x)./x.^2;
I=quadgk(f,2/pi,+Inf)%求得I=1.0000

基于梯形积分法
已知(x_i,y_i)（i=1,2,……,n），且a=x₁＜x₂＜……＜x_n=b，求I=∫abf(x)dx\int_{a}^{b}f(x)dx∫abf(x)dx近似值
I=trapz(x,y)
其中，向量x、y定义函数关系y=f(x)
trapz函数采用梯形积分法则，积分的近似值为：
I=∑i=1n−1hiyi+1+yi2\sum_{i=1}^{n-1}h_i{{y_{i+1}+y_i}\over 2}∑i=1n−1hi2yi+1+yi
其中,h_i=x_i+1-x_i,可以用以下语句实现：
.>>sum(diff(x).*(y(1:end-1)+y(2:end))/2)
例子：设x=1:6，y=[6,8,11,7,5,2]，用trapz函数计算定积分

x=1:6;
y=[6,8,11,7,5,2];
I1=trapz(x,y)%求得I1=35
I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)%求得I2=35

(3)多重定积分的数值求解

求二重积分的数值解：∫cd∫abf(x,y)dxdy\int_{c}^{d}\int_{a}^{b}f(x,y)dxdy∫cd∫abf(x,y)dxdy
I=integral2(filename,a,b,c,d)
I=quad2d(filename,a,b,c,d)
I=dbquad(filename,a,b,c,d,tol)
求三重积分的数值解：∫ef∫cd∫abf(x,y)dxdy\int_{e}^{f}\int_{c}^{d}\int_{a}^{b}f(x,y)dxdy∫ef∫cd∫abf(x,y)dxdy
I=integral3(filename,a,b,c,d,e,f)
I=triplrquad(filename,a,b,c,d,e,f,tol)
例子：∫−11∫−22e−x22sin⁡(x2+y)dxdy\int_{-1}^{1}\int_{-2}^{2}e^{{-x^2\over 2}{\sin(x^2+y)}}dxdy∫−11∫−22e2−x2sin(x2+y)dxdy
∫01∫0π∫0π4xze−z2y−x2dxdydz\int_{0}^{1}\int_{0}^{\pi}\int_{0}^{\pi}4xze^{-z^2y-x^2}dxdydz∫01∫0π∫0π4xze−z2y−x2dxdydz

f1=@(x,y)exp(-x.^2/2).*sin(x.^2+y);
I1=quad2d(f1,-2,2,-1,1)%求得I1=1.5745
f2=@(x,y,z)4*x.*z.*exp(-z.^2.*y-x.^2);
I2=integral3(f2,0,pi,0,pi,0,1)%求得I2=1.7328

二、线性方程求解

1、直接法

高斯（Gauss）消去法
列主元消去法
矩阵的三角分解法

(1)利用左除符号
MATLAB提供了一个左除运算符“\”用于求解线性方程组。如线性方程组Ax=b → x=A^-1b→x=A\b
如果矩阵A是奇异或者接近奇异的，会给出警告信息。
例子：用左除运算符求解下列线性方程组
{2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0\left\{ \begin{aligned} 2x_1+x_2-5x_3+x_4=13\\ x_1-5x_2+7x_4=-9\\ 2x_2+x_3-x_4=6\\ x_1+6x_2-x_3-4x_4=0 \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0

(2)利用矩阵分解求解线性方程组
矩阵分解是将一个给定的矩阵分解成若干个特殊类型矩阵的乘积，从而将一个一般的矩阵计算问题转化为几个易求得特殊矩阵的计算问题。此方法优点是运算速度快，可以节省存储空间。

LU分解
①矩阵的LU分解是将一个n阶矩阵A表示为一个下三角矩阵L和一个上三角矩阵U的乘积。只要方阵是非奇异的，LU分解总可以进行。
Ax=b→Ly=b，Ux=y
②MATLAB提供函数lu，有两种调用格式：
[L,U]=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L，使之满足A=LU。注意，这里的矩阵A必须是方阵。
[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。注意，这里的矩阵A必须是方阵。
③所以，可以利用LU分解求解线性方程组：
Ax=b→LUx=b→x=U(L\b)
或
Ax=b→PAx=Pb→LUx=Pb→x=U(L\P*b)
例子：用LU分解求解下列线性方程组
{2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0\left\{ \begin{aligned} 2x_1+x_2-5x_3+x_4=13\\ x_1-5x_2+7x_4=-9\\ 2x_2+x_3-x_4=6\\ x_1+6x_2-x_3-4x_4=0\\ \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧2x1+x2−5x3+x4=13x1−5x2+7x4=−92x2+x3−x4=6x1+6x2−x3−4x4=0
QR分解
Cholesky分解

2、迭代解法

(1)雅可比（Jacobi）迭代法
对于Ax=b，将A分解成（D-L-U），其中D为对角矩阵，L为负的下三角阵，U为负的上三角阵→(D-L-U)x=b

求解公式为：
x=D−1(L+U)x+D−1bx=D^{-1}(L+U)x+D^{-1}bx=D−1(L+U)x+D−1b
与之对应的迭代公式为：
x(k+1)=D−1(L+U)x(k)+D−1bx^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}bx(k+1)=D−1(L+U)x(k)+D−1b
(令B=D−1(L+U),f=D−1b令B=D^{-1}(L+U),f=D^{-1}b令B=D−1(L+U),f=D−1b)
x(k+1)=Bx(k)+fx^{(k+1)}=Bx^{(k)}+fx(k+1)=Bx(k)+f
雅可比迭代法的函数文件jacobi.m

function [y,n]=jacobi(A,b,x0,ep)%xo为迭代初值,ep为迭代精度
D=diag(diag(A));%生成A的对角阵；
L=-tril(A,-1);%生成A的负的下三角阵
U=-triu(A,1);%生成A的负的上三角阵
B=D\(L+U);
f=D\b;
y=B*x0+f;%根据初值x0,求第一次迭代
n=1;
while norm(y-x0)>=ep%根据前后两次迭代结果的2范数是否很接近来判断x0=y;y=B*x0+f;n=n+1;
end

(2)高斯-赛德尔迭代法(Gauss-Serdel)

迭代公式为：
Dx(k+1)=(L+U)x(k)+bDx^{(k+1)}=(L+U)x^{(k)}+bDx(k+1)=(L+U)x(k)+b
→Dx(k+1)=Lx(k+1)+Ux(k)+b→Dx^{(k+1)}=Lx^{(k+1)}+Ux^{(k)}+b→Dx(k+1)=Lx(k+1)+Ux(k)+b
→(D−L)x(k+1)=Ux(k)+b→(D-L)x{(k+1)}=Ux^{(k)}+b→(D−L)x(k+1)=Ux(k)+b
→x(k+1)=(D−L)−1Ux(k)+(D−L)−1b→x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b→x(k+1)=(D−L)−1Ux(k)+(D−L)−1b
(令B=(D−L)−1U,f=(D−L)−1b令B=(D-L)^{-1}U,f=(D-L)^{-1}b令B=(D−L)−1U,f=(D−L)−1b)
x(k+1)=Bx(k)+fx^{(k+1)}=Bx^{(k)}+fx(k+1)=Bx(k)+f
雅可比迭代法的函数文件GaussSerdel.m

function [y,n]=jacobi(A,b,x0,ep)%xo为迭代初值,ep为迭代精度
D=diag(diag(A));%生成A的对角阵；
L=-tril(A,-1);%生成A的负的下三角阵
U=-triu(A,1);%生成A的负的上三角阵
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;%根据初值x0,求第一次迭代
n=1;
while norm(y-x0)>=ep%根据前后两次迭代结果的2范数是否很接近来判断x0=y;y=B*x0+f;n=n+1;
end

(3)例子

分别用雅可比迭代法和高斯赛德尔迭代法求解方程组。迭代初值为0，精度为10−610^{-6}10−6
[4−2−1−243−1−33][X1X2X3]=[150]\begin{bmatrix} 4 & -2 & -1 \\ -2 & 4 & 3 \\ -1 & -3 & 3 \end{bmatrix} \begin{bmatrix} X_1\\ X_2 \\ X_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 5 \\ 0 \end{bmatrix} ⎣⎡4−2−1−24−3−133⎦⎤⎣⎡X1X2X3⎦⎤=⎣⎡150⎦⎤

A=[4,-2,-1;-2,4,3;-1,-3,3];
b=[1,5,0]';
[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)
[x,n]=GaussSerdel(A,b,[0,0,0]',1.0e-6)

分别用雅可比迭代法和高斯赛德尔迭代法求解方程组。迭代初值为0，精度为10−610^{-6}10−6
[12−2111221][X1X2X3]=[976]\begin{bmatrix} 1 & 2 & -2 \\ 1 & 1 & 1 \\ 2 & 2 & 1 \end{bmatrix} \begin{bmatrix} X_1\\ X_2 \\ X_3 \end{bmatrix} = \begin{bmatrix} 9 \\ 7 \\ 6 \end{bmatrix} ⎣⎡112212−211⎦⎤⎣⎡X1X2X3⎦⎤=⎣⎡976⎦⎤

A=[1,2,-2;1,1,1;2,2,1];
b=[9;7;6];
[x,n]=jacobi(A,b,[0;0;0],1.0e-6
[x,n]=GaussSerdel(A,b,[0;0;0],1.0e-6)

对比上述两个例子，雅可比迭代法和高斯赛德尔迭代法，其是否收敛、收敛速度的快慢，与实际线性方程组的结构有关

3、直接法与迭代法的对比

直接法：以矩阵初等变换为基础，可以求得方程组的精确解；但是占用的空间内存大、程序实现较为复杂；一般适合求解低阶稠密线性方程组。
迭代法：从给定初始值逐步逼近精确值的过程，求解过程占用存储空间小，程序设计简单；适合求解大型稀疏矩阵线性方程组；要考虑算法的收敛性。

三、非线性方程求解与函数极值计算

1、非线性方程数值求解

（1）单变量非线性方程求解
x=fzero(filename,x0)
其中，filename是待求根方程左端的函数表达式，x0是初始值。
①例子：求f(x)=x−1x+5在x0=−5、x0=0.1和x0=1f(x)=x-{1\over x}+5在x_0=-5、x_0=0.1和x_0=1f(x)=x−x1+5在x0=−5、x0=0.1和x0=1作为迭代初值时的根。

f=@(x)x-1./x+5;
x1=fzero(f,-5)
x2=fzero(f,1)
x3=fzero(f,0.1)
fplot(f,[-10,2])
grid on
hold on
fplot(0,[-10,2])
hold on
plot(x1,0,'*',x2,0,'o',x3,0,'p')

求得x1=−5.1926,x2=0.1926,x3=3.7372∗10−16x_1=-5.1926,x_2=0.1926,x_3=3.7372*10^{-16}x1=−5.1926,x2=0.1926,x3=3.7372∗10−16，由图像看，显然x₃不是方程的根。
所以，利用fzero()函数求解方程，初值的选取很重要。
②求f(x)=x2−1f(x)=x^2-1f(x)=x2−1的根

f=@(x)x.^2-1;
x=[];
x0=-0.25:0.001:0.25;
for x00=x0x=[x,fzero(f,x00)];
end
plot(x0,x,'-o');
xlabel('初值');
ylabel('方程的根');
axis([-0.25,0.25,-1,1])

结果表明，相同的根对应的处置范围并不连续，求得的根也并非是离初值比较近的根。所以fzero函数执行的是一个数值搜索的过程，搜索结果依赖于函数特性和指定的函数初值。
（2）非线性方程组的求解
x=fsolve(filename,x0,option)
其中，x为返回的近似解，filename是待求根方程左端的函数表达式，x0是初始值，，option用于设置优化工具箱的优化参数，可以调用optimset函数来完成。例如，Display参数设置为’off’时不显示中间结果。
例子：①利用fsolve函数解上面例子的方程。求f(x)=x−1x+5在x0=−5、x0=0.1和x0=1f(x)=x-{1\over x}+5在x_0=-5、x_0=0.1和x_0=1f(x)=x−x1+5在x0=−5、x0=0.1和x0=1作为迭代初值时的根。

当初值是0.1时，fzero函数无法得到正确的结果，而利用fsolve函数可以。因为不同函数的实现算法不同，适用的场合也不同。
②求下列方程组在(1,1,1)附近的解并对结果进行验证
sinx+y+z2ez=0sinx+y+z^2e^z=0sinx+y+z2ez=0
s+y+z=0s+y+z=0s+y+z=0
xyz=0xyz=0xyz=0
（因为迭代法是逼近问题，为了展示效果，设置为format long）

2、函数极值的计算

函数极值包括极大值和极小值，或者叫最大值和最小值。MATLAB只考虑吧最小值的计算。
（1）无约束最优化问题
minf(x)，其中x=[x1,x2,……,xn]Tminf(x)，其中x=[x_1,x_2,……,x_n]^Tminf(x)，其中x=[x1,x2,……,xn]T
求最小值的函数为：
[xmin,fmin]=fminbnd(filename,x1,x2,option)
[xmin,fmin]=fminsearch(filename,x0,option)
[xmin,fmin]=fminunc(filename,x0,option)
其中，xmin表示极小值点，fmin表示最小值。第一个函数的输入变量x1、x2分别表示被研究区间的左右边界。后两个函数的输入变量x0是一个向量，表示极值点的初值。
例子：求函数f(x)=x−1x+5f(x)=x-{1\over x}+5f(x)=x−x1+5在区间（-10，-1）上的最小值

（2）有约束最优化问题
minf(x)，其中x=[x1,x2,……,xn]T，且使得G(x)≤0minf(x)，其中x=[x_1,x_2,……,x_n]^T，且使得G(x)≤0minf(x)，其中x=[x1,x2,……,xn]T，且使得G(x)≤0
即求取一组x，使得目标函数f(x)为最小，且满足约束条件G(x)≤0。
函数为：
[xmin,fmin]=fmincon(filename,x0,A,b,Aeq,beq,Lb,Ub,nonlcon,option)
其中，A,b,Aeq,beq,Lb,Ub,nonlcon分别表示线性不等式约束、线性等式约束、x的下界和上界、定义非线性约束的函数。如果某个约束不存在，用空矩阵表示。
具体各种约束，请查看：https://ww2.mathworks.cn/help/optim/ug/fmincon.html
例子：①求解有约束的最优化问题
minf(x)=0.4x2+x12+x22−x1x2+130x12，其中x1,x2,x3满足以下约束条件：minf(x)=0.4x_2+x_1^2+x_2^2-x_1x_2+{1\over 30}x_1^2，其中x_1,x_2,x_3满足以下约束条件：minf(x)=0.4x2+x12+x22−x1x2+301x12，其中x1,x2,x3满足以下约束条件：
x1+0.5x2≥0.4，0.5x1+x2≥0.5，x1≥0，x0≥0x_1+0.5x_2≥0.4，0.5x_1+x_2≥0.5，x_1≥0，x_0≥0x1+0.5x2≥0.4，0.5x1+x2≥0.5，x1≥0，x0≥0

f=@(x)0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;
x0=[0.5;0.5];
A=[-1,-0.5;-0.5,-1];%Ax≤b
b=[-0.4;-0.5];
lb=[0;0];%lb≤x≤ub
option=optimset('Display','off');
[xmin,fmin]=fmincon(f,x0,A,b,[],[],lb,[],[],option)

例子：②此例特殊说明非线性约束
（x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) 执行最小化时，满足 nonlcon 所定义的非线性不等式 c(x) 或等式 ceq(x)。fmincon 进行优化，以满足 c(x) ≤ 0 和 ceq(x) = 0。）

编写函数non.m定义非线性约束条件:

function [c,ceq]=non(x);
c=[-x(1).^2+x(2)-x(3).^2x(1)+x(2).^2+x(3).^3-20];
ceq=[-x(1)-x(2).^2+2x(2)+2*x(3).^2-3];

编写主程序函数

f=@(x)x(1).^2+x(2).^2+x(3).^2+8;;
[xmin,fmin]=fmincon(f,rand(3,1),[],[],[],[],zeros(3,1),[],'non')

（3）最小值问题实例

①假设仓库所选点坐标为（x，y），则总里程表达式为：
d(x,y)=10(x−10)2+(y−10)2+18(x−30)2+(y−50)2+20(x−16.667)2+(y−29)2+14(x−0.555)2+(y−29.888)2+25(x−22.2221)2+(y−49.988)2d(x,y)=10\sqrt{(x-10)^2+(y-10)^2}+18\sqrt{(x-30)^2+(y-50)^2}+20\sqrt{(x-16.667)^2+(y-29)^2}+14\sqrt{(x-0.555)^2+(y-29.888)^2}+25\sqrt{(x-22.2221)^2+(y-49.988)^2}d(x,y)=10(x−10)2+(y−10)2+18(x−30)2+(y−50)2+20(x−16.667)2+(y−29)2+14(x−0.555)2+(y−29.888)2+25(x−22.2221)2+(y−49.988)2

a=[10,30,16.667,0.555,22.2221];
b=[10,50,29,29.888,49.988];
c=[10,18,20,14,25];
f=@(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
[xmin,fmin]=fminsearch(f,[15,30])

②若地域限制，仓库必须建立在曲线y=x2y=x^2y=x2上

编写函数non.m定义非线性约束条件:

function [c,ceq]=non(x);
c=[];
ceq=[x(2)-x(1)^2];

编写主程序函数

a=[10,30,16.667,0.555,22.2221];
b=[10,50,29,29.888,49.988];
c=[10,18,20,14,25];
f=@(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
option=optimset('Display','off');
[xmin,fmin]=fmincon(f,[15,30],[],[],[],[],[],[],'non',option)

四、常微分方程数值求解

1、常微分方程数值求解的一般概念

求解常微分方程初值问题就是寻找函数y(t)使之满足如下方程：
{y′=f(t,y)，t0≤t≤by(t0)=y0\left\{ \begin{aligned} y'=f(t,y)，t_0≤t≤b\\ y(t_0)=y_0 \end{aligned} \right. {y′=f(t,y)，t0≤t≤by(t0)=y0
所谓其数值解法，就是求y(t)在离散节点t_n处的函数近似值yny_nyn的方法，yn≈y(xn)y_n≈y(x_n)yn≈y(xn)。相邻两个节点之间的距离称为步长。一般有单步法和多步法。

2、常微分方程数值求解函数

（1） [t,y]=solver(filename,tspan,y0,option)
其中，t和y分别给出时间向量和相应的数值解。solver为求常微分方程数值解的函数。filename为定义f(t,y)的函数名，该函数必须返回一个列向量。tspan形式为[t0,tf]，表示求解区间。y0是初始状态向量。option是可选参数，用于设置求解属性。常用的属性包括相对误差值RelTol（默认值是10−310^{-3}10−3）和绝对误差值AbsTol（默认值是10−610^{-6}10−6）。
（2）常微分方程数值求解函数的统一命名格式：
odennxx

ode是Ordinary Differential Equation
nn是数值，代表使用方法的阶数
xx是字母，用于标注方法的专门特征

例子：求微分方程初值问题，并与精确解y(t)=t+1+1y(t)=\sqrt{t+1}+1y(t)=t+1+1进行比较
{y′=y2−t−24(t+1)y(0)=2，0≤t≤10\left\{ \begin{aligned} y'={y^2-t-2}\over {4(t+1)} \\ y(0)=2 \end{aligned} \right. ，0≤t≤10⎩⎪⎨⎪⎧4(t+1)y′=y2−t−2y(0)=2，0≤t≤10

f=@(t,y)(y^2-t-2)/4/(t+1);
[t,y]=ode23(f,[0,10],2);
y1=sqrt(t+1)+1;
plot(t,y,'b:',t,y1,'r');

例子：已知一个二阶线性系统的微分方程为：
{d2xdt2+ax=0，a>0x(0)=0，x′(0)=1\left\{ \begin{aligned} {{d^2}x \over {dt^2}}+ax=0 ，a>0\\ x(0)=0，x'(0)=1 \end{aligned} \right. ⎩⎨⎧dt2d2x+ax=0，a>0x(0)=0，x′(0)=1
取a=2，绘制系统的时间相应曲线和相平面图。
对于高阶的微分方程，一般需要进行降阶处理。令x2=x，x1=x′x_2=x，x_1=x'x2=x，x1=x′，则得到系统的状态方程：
{x1′=−ax2x2′=x1x2(0)=0,x1(0)=1\left\{ \begin{aligned} x_1'=-ax_2\\ x_2'=x_1\\ x_2(0)=0,x_1(0)=1 \end{aligned} \right. ⎩⎪⎨⎪⎧x1′=−ax2x2′=x1x2(0)=0,x1(0)=1

f=@(t,x)[0,-2;1,0]*[x(1);x(2)];
[t,x]=ode45(f,[0,20],[1,0]);
subplot(1,2,1);
plot(t,x(:,2));%x(:,2)就是x(2),x(2)=x
subplot(1,2,2);
plot(x(:,2),x(:,1));

3、刚性问题

刚性问题是指一类微分方程，其解的分量有的变化很快，有的变化很慢，且相差悬殊。
刚性问题的数值解算法必须取很小步长才能获得满意的结果。解决刚性问题需要专门的方法，虽然非刚性算法也可以求解，但是需要很长的计算时间。

例子：

当lamda=0.01时

lamda=0.01;
f=@(t,y)y^2-y^3;
tic;%启动计时器
[t,y]=ode45(f,[0,2/lamda],lamda);
toc;%结束计时器
disp(['ode45计算的点数' num2str(length(t))]);

历时 0.017177 秒。
ode45计算的点数157

当lamda=1e-5时

lamda=1e-5;
f=@(t,y)y^2-y^3;
tic;%启动计时器
[t,y]=ode45(f,[0,2/lamda],lamda);
toc;%结束计时器
disp(['ode45计算的点数' num2str(length(t))]);

历时 0.392539 秒。
ode45计算的点数120565
这时计算时间明显加长，计算的点数剧增，表明此时常微分方程表现为刚性

对于刚性问题，我们选择以“s”结尾的函数

lamda=1e-5;
f=@(t,y)y^2-y^3;
tic;%启动计时器
[t,y]=ode15s(f,[0,2/lamda],lamda);
toc;%结束计时器
disp(['ode15s计算的点数' num2str(length(t))]);

历时 0.770228 秒。
ode15s计算的点数98
很明显，计算时间大大缩短，计算的点数大大减少。

五、常微分方程应用举例

1、Lotka-Volterra模型

（1）模型分析
x1(t)x_1(t)x1(t)—兔子在t时刻的数量，x2(t)x_2(t)x2(t)—狐狸在t时刻的数量
r1r_1r1—兔子独自生存时的增长率，r2r_2r2—狐狸独自存在时的死亡率
λ1(t)λ_1(t)λ1(t)—狐狸掠取兔子的能力系数，λ2(t)λ_2(t)λ2(t)兔子对狐狸的供养能力系数
若假设r1=2，r2=1，λ1=λ2=0.01r_1=2，r_2=1，λ_1=λ_2=0.01r1=2，r2=1，λ1=λ2=0.01则，兔子的增长模型如下：
dx1dt=x1(r1−λ1x2){dx_1\over dt}=x_1(r_1-λ_1x_2)dtdx1=x1(r1−λ1x2)
狐狸的增长模型如下：
dx2dt=x2(−r1+λ2x1){dx_2\over dt}=x_2(-r_1+λ_2x_1)dtdx2=x2(−r1+λ2x1)
（2）第一问：兔子数量初始值x1(0)=300x_1(0)=300x1(0)=300，狐狸数量初始值x2(0)=150x_2(0)=150x2(0)=150。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150]);
subplot(1,2,1);
plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('x1(t)','x2(t)');
xlabel('时间');ylabel('物种数量');
grid on
subplot(1,2,2);
plot(x(:,1),x(:,2));
grid on

结论：兔子和狐狸两者相互制约，当狐狸逐渐增多时，兔子逐渐减少。当狐狸数量增加到一定数量时，由于种群内部竞争激烈，又导致其数量减少，从而兔子天敌减少，导致兔子数量增加。当兔子数量增加到一定数量时，由于种群内部竞争激烈，又导致其数量减少。如此循环，相互制约发展。
（3）第二问：兔子数量初始值x1(0)=15x_1(0)=15x1(0)=15，狐狸数量初始值x2(0)=22x_2(0)=22x2(0)=22。

（4）第三问：兔子数量初始值x1(0)=102x_1(0)=102x1(0)=102，狐狸数量初始值x2(0)=198x_2(0)=198x2(0)=198。

（5）验证（1/λ，2/λ）是稳定平衡点

取λ=0.01，所以稳定平衡点（1/λ，2/λ）即（100，200），以此点作为初值画图。
进一步分析初始值变化之后的图像
①当初始值变为(70,150)，向下偏离平衡点比较远，其图像：

可发现兔子和狐狸的数量变化比较剧烈，但是依然成周期性变化。
②当初始值变为(900,1600)，向上偏离平衡点很远，其图像：

可发现，兔子和狐狸的数量变化十分剧烈，但是依然成周期性变化。

2、Lotka-Volterra改进模型

（考虑环境的承受力，限制兔子的数量上限）

（1）第一问：原模型下，狐狸和兔子数量的时间函数曲线

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');ylabel('物种数量');
title('原模型下，狐狸和兔子数量的时间函数曲线')
grid on

（2）第二问：改进模型下，狐狸和兔子数量的时间函数曲线（令R=400，λ=0.01）

rabbitFox=@(t,x)[2*x(1)*(1-x(1)/400-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');ylabel('物种数量');
title('改进模型下，狐狸和兔子数量的时间函数曲线')
grid on

原模型无论经历多长时间，狐狸和兔子的数量总是在自己的平衡点之间波动。而改进之后的模型，虽然在前一段时间有较大的波动，但随着时间的推移，波动越来越小。在经历足够长的时间后，狐狸和兔子的数量分别达到稳定平衡。这更加接近自然界的实际情况。

（3）第三问：原模型下，狐狸数量相对于兔子数量的函数曲线。

rabbitFox=@(t,x)[x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150]);
plot(x(:,1),x(:,2));
xlabel('时间');ylabel('物种数量');
grid on
title('原模型下，狐狸数量相对于兔子数量的函数曲线')

（4）第四问：改进模型下，狐狸数量相对于兔子数量的函数曲线。

rabbitFox=@(t,x)[2*x(1)*(1-x(1)/400-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(x(:,1),x(:,2));
grid on
xlabel('时间');ylabel('狐狸数量');
title('改进模型下，狐狸数量相对于兔子数量的函数曲线。')
grid on

总结

专题六数值微积分与方程求解相关推荐

数学建模之matlab软件学习06——专题六数值微积分与方程求解
#本文章仅用于记录本人学习过程,当作笔记来用,如有侵权请及时告知,谢谢! 6.1 数值微分与数值积分数值微分 6.2 6.3 线性方程组应用举例: 平面桁架结构受力分析问题小行星运行轨道计算问题: ...
6.数值微积分与方程求解
数值微分与数值积分 % 数值积分 (有函数式版本定积分) 方法1 % [I,n]=quad(filename,a,b,tol,trace) 基于自适应辛普森方法 % [I,n]=quadl(file ...
数值微积分与方程求解
1.数值微分与数值积分数值微分 MATLAB提供了求向前差分的函数diff,其调用格式有三种: dx=diff(x):计算向量x的一阶向前差分,dx(i)=x(i+1)-x(i),i=1,2,-,n ...
6. 数值微积分与方程求解
文章目录 A 数值微分与数值积分 A.a数值微分(diff) A.b 数值积分 B 线性方程组求解 B.a 直接法 B.b 迭代法 B.c 线性方程组应用举例 C 非线性方程求解与函数极值计算 C.a ...
matlab积分练习,matlab练习之数值微积分和方程数值求解
一.符号积分求符号积分函数:int 格式:int(f,x,a,b) 功能:计算定积分格式:int(f,x) 功能:计算不定积分使用int函数之前,先用syms声明x是符号变量例: 代码: sy ...
非线性规划MATLAB求解原理,专题六--非线性规划介绍及其Matlab求解方法.ppt
迭代法一般步骤注意:数值求解最优化问题的计算效率取决于确定搜索方向P (k)和步长的效率. Matlab求解方法简介 Step3: 利用(3)式或其它一维搜索的方法求计算然后令k:=k+1, ...
工程数学（数值分析）第六讲：数值微积分
文章目录第六讲:数值微积分二点.三点数值微积分求导数定步长梯形求积公式变步长梯形求积公式定步长.变步长辛普森求积公式第六讲:数值微积分二点.三点数值微积分求导数定步长梯形求积公式变步 ...
matlab整理符号表达式,[2018年最新整理]MATLAB符号运算与符号方程求解.ppt
[2018年最新整理]MATLAB符号运算与符号方程求解 MATLAB符号计算 1 符号对象 2 符号微积分 3 级数 4 符号方程求解 9.1 符号对象 9.1.1 建立符号对象 1．建立符号变量 ...
偏微分方程数值解法python_Python数值计算----------求解简单的偏微分方程
很多物理现象的都可以用方程来描述,比如热传导与物质扩散可以用扩散方程来描述,流体的流动可以用NS方程描述等等.如果能够将这些偏微分方程求解出来,就可以来对很多物理现象进行仿真,现在工程中的仿真软件都是 ...

专题六数值微积分与方程求解