如何利用matlab求解方程
1.    前言
作为三大数学软件之一，matlab在数值计算方法的能力首屈一指。求解方程是工科学习和工程计算中最基础、最常见的问题。掌握利用现代化工具求解方程的方法对于提升我们的工科素养至关重要。为此，本篇将对matlab中求解方程的方法进行介绍。
2.    用法
求解过程
2.1    指明变量
告诉电脑方程中所含有的变量，包括参数和未知变量。比如：所求解的方程为：，很显然该方程中有a,b,c,x符号变量，因此该步骤的写法为：
  
syms a b c x  
  
 
 
2.2    指明方程，未知数和限制条件(非必需)
eqns
  
方程，如果超过一个，则放在[ ]中，并用逗号隔开。如：
  
 
Vars
  
待求的未知数
 
Names-value(非必需)
  
Names:‘returnConditions’
  
是否返回出含有参数的通解。’true‘为返回，’false‘为否，即给出一个特解；
 
Name: 'IgnoreAnalyticConstraints'
  
是否给出解的最简形式。
  
‘true‘为是，‘false’为否
 
Name:'PrincipalValue'
  
是否仅给出一个解。False为返回所有的解，true为仅返回一个解；
 
Name:’Real’
  
是否仅返回实数解
 
2.3    获得所求方程的解
如果为多个函数，该解存储形式为结构体。
 
3.    具体实例
3.1   求解sin(x)=1的通解
 
具体代码：
syms x  
  
[x,params,conds]=solve(sin(x)==1,'ReturnConditions', true)  
 
 
结果
solx =pi/2+2*pi*k   
  
params =k   
  
conds =in(k,'integer')
 
 
可以看出，该方程的通解为：
 
3.2   求解以下方程：
代码：
syms a b c y x
  
[x,y]=solve([a*x^2+b*y+c==0,a*x+2*y==4],[x,y])
 
 
结果：
x =
  
 ((a*b)/4-(-(a*(- a*b^2+32*b +16*c))/16)^(1/2))/a
  
 ((a*b)/4+(-(a*(- a*b^2+32*b +16*c))/16)^(1/2))/a
  
 
  
y =
  
 (-(a*(- a*b^2+32*b +16*c))/16)^(1/2)/2-(a*b)/8+2
  
 2-(-(a*(- a*b^2+32*b +16*c))/16)^(1/2)/2-(a*b)/8
 
 
即：

文章目录
前言
1 直接法
2 迭代法
小结

前言
今天我们要说的就是数值微积分，赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第2小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。
1 直接法
1.线性方程组的直接解法

高斯（Gauss）消去法

列主元消去法

矩阵的三角分解法

高斯（Gauss）消去法是一个经典的直接法，由它改进得到的列主元消去法，是目前计算机上求解线性方程组的标准算法，其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题。此外，还有矩阵的三角分解法等许多直接求解算法。

（1）利用左除运算符的直接解法

MATLAB提供了一个左除运算符“\”用于求解线性方程组，它使用列主元消去法，使用起来十分方便。对于线性方程组$Ax=b$，可以利用左除运算符反斜杠求解，b左除以A可获得线性方程组的数值解x。

$Ax=b$→$x=A\backslash b$

注意，如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。

例1  用左除运算符求解下列线性方程组。

$\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.$
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
x=A\b

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

LU分解

QR分解

Cholesky分解

①LU分解的基本思想

矩阵的LU分解就是将一个n阶矩阵表示为一个下三角矩阵和一个上三角矩阵的乘积。线性代数中已经证明，只要方阵是非奇异的，LU分解总是可以进行的。

$A=LU →Ax=b→LUx=b$→ $\left\{
\begin{aligned}Ly=b\\Ux=y
\end{aligned}
\right.$

对于三角方程很容易求解，于是可以首解向量y使Ly=b，再求解Ux=y，从而达到求解线性方程组Ax=b的目的。

②MATLAB的LU分解函数

LU分解函数是根据列主元LU分解算法定义的，具有较好的数据稳定性。lu函数有两种调用格式：

[L，U]=lu（A）：

产生一个上三角阵U和一个变换形式的下三角阵L，使之满足A=LU。注意，这里的矩阵A必须是方阵。

[L，U，P]=lu（A）：

产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。同样，矩阵A必须是方阵。

当使用第一种格式时，矩阵L往往不是一个下三角阵，但可以通过行交换成为一个下三角阵。

③用LU分解求解线性方程组

$Ax=b$→$LUx=b$→$x=U \backslash (L\backslash b）$

或

$Ax=b$→$PAx=Pb$→$LUx=Pb$→$x=U \backslash（L*b）$

通过LU分解后可以大大提高运算速度。

例2  用LU分解求解例1中的线性方程组。
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U\(L\b)

x =
  -66.5556
   25.6667
  -18.7778
   26.5556

2 迭代法
2.线性方程组的迭代解法

迭代法是一种不断用变量的原值推出它的新值的过程，是用计算机解决问题的一种基本方法。

$\left\{
\begin{aligned}
10x_1-x_2=9\\
-x_1+10x_2-2x_3=7\\
-2x_2+10x_3=6
\end{aligned}
\right.$$\implies$ $\left\{
\begin{aligned}
x_2&=10x_1-9\\
x_1&=10x_2-2x_3-7\\
x_3&=(6+2x_2)10
\end{aligned}
\right.$$\implies$$\left\{
\begin{aligned}
x_1^{k+1}&=10x_2^{k}-2x_3^{k}-7\\
x_2^{k+1}&=10x_1^{k}-9\\
x_3^{k+1}&=0.6+0.2x_2^{k}
\end{aligned}
\right.$

（1）雅可比（Jacobi）迭代法

$Ax=b,
A=D-L-U，
(D-L-U)x=b$

$D=\left [
\begin{array}{cccc}
 a_{11} & & & \\
               & a_{22} &  & \\
                &  & \ddots & \\
                                          & &   & a_{nn} 
\end{array}
\right]$

$L=-\left [
\begin{array}{cccc}
  0 &  &&  &  \\
  a_{21}  & 0 &&  & \\
  a_{31}  &  a_{32}   &0&  &  \\
\vdots &\vdots&\vdots &\ddots&\\
  a_{n1} &  a_{n2}&\cdots& \cdots&0
\end{array}
\right]$

$U=-\left [
\begin{array}{cccc}
  0 & a_{12} &a_{13}& \cdots & a_{1n} \\
               & 0 &a_{23}& \cdots &a_{2n} \\
                &  &\ddots& \ddots & \vdots \\
                          &   & &\ddots&a_{n-1n}\\
                                          &&&   &0
\end{array}
\right]$
求解公式为：

$x=D^{-1}（L+U）x+D^{-1}b$

与之对应的迭代公式为：

$x^{k+1}=D^{-1}（L+U）x^{k}+D^{-1}b$ $\implies(B=D^{-1}(L+U),f=D^{-1}b$) $\implies x^{k+1}=Bx^{k}+f$

雅可比迭代法的函数文件jacobi.m：
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
    x0=y;
    y=B*x0+f;
    n=n+1;
end

（2）高斯-赛德尔（Gauss-Serdel）迭代法

$Dx^{k+1}=（L+U）x^{k}+b
→Dx^{k+1}=Lx^{k+1}+Ux^k+b→
（D-L）x^{k+1}=Ux^{k}+b \\
x^{k+1}=（D-L）^{-1}U×0+（U-L）^{-1}b
→B=D-L^{-1}U,f=D-L^{-1}b→ x^{k+1}=Bx^{k}+f$

Gauss-Serdel迭代法的函数文件gauseidel.m
function [y,n]=gauseidel(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1); 
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
    x0=y;
    y=B*x0+f;
    n=n+1;
end



例3  分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0，迭代精度为10^(-6)。

$\left[\begin{matrix}
4&-2&-1\\
-2&4&3\\
-1&-3&3
\end{matrix}\right]$$\left[\begin{matrix}
x_1\\
x_2\\
x_3
\end{matrix}\right]$=$\left[\begin{matrix}
1\\
5\\
0
\end{matrix}\right]$
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]=gauseidel(A,b,[0,0,0]',1.0e-6)

x =
    0.9706
    0.8529
    1.1765
n =
    35
x =
    0.9706
    0.8529
    1.1765
n =
    16

例4  分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组，看是否收敛。

$\left[\begin{matrix}
1&-2&2\\
1&1&1\\
2&2&1
\end{matrix}\right]$$\left[\begin{matrix}
x_1\\
x_2\\
x_3
\end{matrix}\right]$=$\left[\begin{matrix}
9\\
7\\
6
\end{matrix}\right]$
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]=gauseidel(A,b,[0;0;0],1.0e-6)

x =
   -27
    26
     8
n =
     4
x =
   NaN
   NaN
   NaN
n =
    1012

小结
直接法：以矩阵初等变换为基础，可以求得方程组的精确解；占用的内存空间大、程序实现较为复杂；一般适合求解低阶稠密线性方程组。

迭代法：从给定初始值逐步逼近精确解的过程，求解过程占用存储空间小，程序设计简单；适用于求解大型稀疏矩阵线性方程组；要考虑算法的收敛性。

大家学会了吗？另外分享给大家一个Markdown的公式神器-$\ Katex$ ，为了编辑出优美的公式，大家可以多看看官方支持文档.最后没有一键三连，但欢迎大家

点赞👍，收藏⭐，转发🚀，

如有问题、建议，请您在评论区留言💬哦。

求微分方程的解
自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。下面介绍如何用 Matlab 来计算微分方程（组）的数值解。
Euler折线法
微分方程的基本数值解法——Euler折线法。
步骤：

分割求解区间，差商代替微商，解代数方程
例子：

解：
等距剖分：$a=x_0<x_1<x_2<	\cdots <x_{n-1}<x_n=b$
步长：$h=x_{k+1}-x_k=(b-a)/n,k=0,1,2,...,n-1$
$y(x_{k+1})\approx y(x_k)+hy\prime(x_k)$
取步长 h = (2 - 0)/n = 2/n，得差分方程



当 h=0.4，即 n=5 时，代码如下：
clear
a=0;  b=2;
h=0.4;
n=(b-a)/h+1; 
x(1)=0;  y(1)=1;
for i=1:n-1        
    y(i+1)=y(i)+h*(y(i)+2*x(i)/y(i)^2);
    x(i+1)=x(i)+h;
end
y1=((5/3)*exp(3*x)-2*x-2/3).^(1/3);
plot(x,y1,x,y);
h = legend('解析解','Euler解');


解得：

$y=(\frac{5}{3}e^{3x}-2x-\frac{2}{3})^{1/3}$
Runge-Kutta方法
龙格－库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
同样是上面的例子:




clear;
a=0; b=2; h=0.4;
n=(b-a)/h+1; 
x(1)=0; y(1)=1;  
for i=1:n-1
    L1=y(i)+2*x(i)/y(i)^2;
    t=x(i)+h/2;t1=y(i)+h*L1/2;
    L2=t1+2*t/t1^2;t1=y(i)+h*L2/2;
    L3=t1+2*t/t1^2;t1=y(i)+h*L3/2;
    L4=t1+2*(x(i)+h)/t1^2;
    y(i+1)=y(i)+h*(L1+2*L2+2*L3+L4)/6; 
    x(i+1)=x(i)+h;
end
plot(x,y)


MATLAB自带的ODE求解器
ODE （常微分方程）
没有一种算法可以有效地解决所有的 ODE 问题，因此MATLAB 提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。



基本语法

[T,Y] = solver(odefun,tspan,y0)

其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T (向量) 中返回的是分割点的值(自变量)，Y (向量) 中返回的是解函数在这些分割点上的函数值。

solver 为Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）

文章目录
前言
1 常微分方程
1.1 常微分方程的概念
1.2 常微分方程的定义
1.3 常微分方程数值求解的一般概念
2 常微分方程数值求解函数
3 刚性问题
结语

前言
今天我们要说的就是利用MATLAB对常微分方程进行数值求解。本文是科学计算与MATLAB语言专题六第5小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。
1 常微分方程
1.1 常微分方程的概念
在初等数学中就有各种各样的方程，比如线性方程、二次方程、高次方程、指数方程、对数方程、三角方程和方程组等等。这些方程都是要把研究的问题中的已知数和未知数之间的关系找出来，列出包含一个未知数或几个未知数的一个或者多个方程式，然后取求方程的解。

但是在实际工作中，常常出现一些特点和以上方程完全不同的问题。比如：物质在一定条件下的运动变化，要寻求它的运动、变化的规律；某个物体在重力作用下自由下落，要寻求下落距离随时间变化的规律；火箭在发动机推动下在空间飞行，要寻求它飞行的轨道，等等，要以现有数据求得出形式上的函数解析式，而不是以已知函数来计算特定的未知数。

物质运动和它的变化规律在数学上是用函数关系来描述的，因此，这类问题就是要去寻求满足某些条件的一个或者几个未知函数。也就是说，凡是这类问题都不是简单地去求一个或者几个固定不变的数值，而是要求一个或者几个未知的函数。

解这类问题的基本思想和初等数学解方程的基本思想很相似，也是要把研究的问题中已知函数和未知函数之间的关系找出来，从列出的包含未知函数的一个或几个方程中去求得未知函数的表达式。但是无论在方程的形式、求解的具体方法、求出解的性质等方面，都和初等数学中的解方程有许多不同的地方。

在数学上，解这类方程，要用到微分和导数的知识。因此，凡是表示未知函数的导数以及自变量之间的关系的方程，就叫做微分方程。
1.2 常微分方程的定义
定义1：凡含有参数，未知函数和未知函数导数 (或微分) 的方程，称为微分方程，有时简称为方程，未知函数是一元函数的微分方程称作常微分方程，未知函数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数，称为微分方程的阶。定义式如下：

$F(x,y,y',y'',…,y^{(n)})=0$ 　

定义2：任何代入微分方程后使其成为恒等式的函数，都叫做该方程的解。

若微分方程的解中含有任意常数的个数与方程的阶数相同，且任意常数之间不能合并，则称此解为该方程的通解(或一般解).

当通解中的各任意常数都取特定值时所得到的解，称为方程的特解。

一般地说，$n$阶微分方程的解含有 $n$个任意常数。

也就是说，微分方程的解中含有任意常数的个数和方程的阶数相同，这种解叫做微分方程的通解。通解构成一个函数族。

如果根据实际问题要求出其中满足某种指定条件的解来，那么求这种解的问题叫做定解问题，对于一个常微分方程的满足定解条件的解叫做特解。

对于高阶微分方程可以引入新的未知函数，把它化为多个一阶微分方程组。
1.3 常微分方程数值求解的一般概念
求解常微分方程初值问题就是寻找函数$y(t)$使之满足如下方程：

$y'=f(t，y),t_0≤t≤b, y(t_0)=y_0$

所谓其数值解法，就是求$y(t)$在离散结点t处的函数近似值$y$。的方法，$y_n\approx y_{xn}$。

这些近似值称为常微分方程初值问题的数值解。相邻两个结点之间的距离称为步长。

单步法：在计算$y_{n+1}$时只用到前一步的$y_n$，因此在有了初值之后就可以逐步往下计算，其代表是龙格-库塔（Runge-Kutta）法。$$

多步法：在计算$y_{n+1}$时，除了用到前一步的值$y_n$之外，还要用到$y_{n-p}（p=1,2,…,k,k>0）$的值，即前面的k步。其代表就是亚当斯$（Adams）$法。
2 常微分方程数值求解函数
MATLAB提供了多个求常微分方程初值问题数值解的函数，一般调用格式为：

[t，y]=solver（filename，tspan，y0，option）

其中，t和y分别给出时间向量和相应的数值解。

solver为求常微分方程数值解的函数。

filename是定义f（t，y）的函数名，该函数必须返回一个列向量。

tspan形式为[to，tf]，表示求解区间。

y0是初始状态向量。

option是可选参数，用于设置求解属性，常用的属性包括相对误差值Re1To1（默认值是$10^{-3}$）和绝对误差值AbsTol（默认值是$10^{-6}$）。

常微分方程数值求解函数的统一命名格式：odennxx

ode是0rdinary Differential Equation的缩写，是常微分方程的意思。

nn是数字，代表所用方法的阶数。例如，ode23采用2阶龙格-库塔$（Runge-Kutta）$算法，用3阶公式做误差估计来调节步长，具有低等精度。

ode45采用4阶龙格-库塔算法，用5阶公式做误差估计来调节步长，具有中等精度。

xx是字母，用于标注方法的专门特征。

ode15s、ode23s中的“s”代表（Stiff），表示函数适用于刚性方程。

求常微分方程数值解的函数




求解函数
采用方法
适用场合




ode23
2阶或3阶Runge-Kutta算法，低精度
非刚性


ode45
4阶或5阶Runge-Kutta算法，中精度
非刚性


ode113
Adams算法，精度可到10-3~10-6
非刚性，计算时间比ode45短


ode15s
Gear’s反向数值微分算法，中精度
刚性


ode23s
2阶Rosebrock算法，低精度
刚性，当精度较低时，计算时间比ode15s短


ode23tb
梯形算法，低精度
刚性，当精度较低时，计算时间比ode15s短


例1求解微分方程初值问题，并与精确解$y（t）=\sqrt {t+1}+1$进行比较。

$\left\{
\begin{aligned}
y'&=\frac{y^2-t-2}{4(t+1)},0≤t≤10\\
y(0)&=2
\end{aligned}
\right.$
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');



例2已知一个二阶线性系统的微分方程为：

$\left\{
\begin{aligned}
\frac{d^2x}{dt^2} +ax&=0,a>0\\
x(0)&=0
,x'(0)=1
\end{aligned}
\right.$

取a=2，绘制系统的时间响应曲线和相平面图。令x2=x，x1=x’，则得到系统的状态方程：

$\left\{
\begin{aligned}
x_2'&=x_1\\
x_1'&=-ax_2\\
x_2(0)&=0,x_1(0)=1\\
\end{aligned}
\right.$
3 刚性问题
有一类常微分方程，其解的分量有的变化很快，有的变化很慢，且相差悬殊，这就是所谓的刚性问题（Stiff）。

对于刚性问题，数值解算法必须取很小步长才能获得满意的结果，导致计算量会大大增加。解决刚性问题需要有专门方法。非刚性算法可以求解刚性问题，只不过需要很长的计算时间。
f=@(t,x)[-2,0;0,1]*[x(2);x(1)];
[t,x]=ode45(f,[0,20],[1,0]); 
subplot(1,2,1);
plot(t,x(:,2)); 
subplot(1,2,2);
plot(x(:,2),x(:,1));



例3假如点燃一个火柴，火焰球迅速增大直至某个临界体积，然后，维持这一体积不变，原因是火焰球内部燃烧耗费的氧气和从球表面所获氧气达到平衡。其简化模型如下：

$\left\{
\begin{aligned}
y'&=y2-y3\\
y(0)&=\lambda 
,0≤t≤2/\lambda
\end{aligned}
\right.$

其中，$y（t）$代表火焰球半径，初始半径是$\lambda$，它很小。分析入的大小对方程求解过程的影响。
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.005774秒。
ode45计算的点数157

tic和toc函数用来记录微分方程求解命令执行的时间，使用tic函数启动计时器，使用toc函数显示从计时器启动到当前所经历的时间。最后还输出计算的点数，运行结果表明这时常微分方程不算很刚性。

 注意：

1 这里的代码易出现错误。例如1和l。

2 $\lambda$的拼写应为lambda。


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))]);

时间已过 1.729913 秒。
ode45计算的点数120565

这时计算时间明显加长，计算的点数剧增，表明这时常微分方程表现为刚性。
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.321166 秒。

ode15s计算的点数98

对于刚性问题，我们需要改变求解算法。我们选择以“s”结尾的函数，他们专门用于求解刚性方程。计算时间大大缩短，计算的点数大大减少。
结语
欢迎大家点赞👍，收藏⭐，转发🚀，

如有问题、建议，请您在评论区留言💬哦。



百度百科-常微分方程。 ↩︎

百度百科-常微分方程。 ↩︎