28、MATLAB线性代数与常微分方程问题解析

MATLAB线性代数与常微分方程问题解析

1. 线性代数练习

1.1 矩阵操作练习

• 获取5x5单位矩阵 ：从5x5的幻方矩阵中获取5x5的单位矩阵。
• 创建10x10单位矩阵 ：从10x10的随机方阵中创建10x10的单位矩阵。
• 复制矩阵 ：从3x3的pascal矩阵中获取3x9的复制方阵。

1.2 方程求解练习

求解以下方程组，并讨论a = 13和a = 29两种情况的解：

q1 + q2 = a
q1 - q2 = 1
2q1 + 3q2 = 23

编写一个带有逻辑和循环运算符（if、break、for和end）的脚本，找出使这些方程有实数解的a值，a为1到50之间的整数。提示：使用rank()函数和反斜杠()运算符。

1.3 多项式求解练习

使用roots()、solve()、zero()和Simulink模型求解以下多项式：
- (2u^9 - 3u^8 + 5u^4 - 13u^2 - 131 = 0)
- (y^7 + 5y^5 - 4y^4 - 13y^3 + 11y^2 + 9y + 3 = 0)
- (5x^4 + 4x^3 + 11x^2 - 3x - 13 = 269)

1.4 数组创建练习

创建一个从10开始到100结束的对数间隔数组（行向量）B，并从行向量B创建列向量BB。

1.5 声音播放练习

播放由以下表达式定义的声音：
[S(t) = \cos(2\pi f_1t) + \cos(2\pi f_2t) + \cot(2\pi f_3t) + \tan(2\pi f_4t) + \tan(2\pi f_5t)]
其中，(f_s = 10000 Hz)（采样频率）；(t = 13 s)（时间长度）；(f_1 = 100 Hz)（第1个信号）；(f_2 = 200 Hz)（第2个信号）；(f_3 = 300 Hz)（第3个信号）；(f_4 = 600 Hz)（第4个信号）；(f_5 = 700 Hz)（第5个信号）。

1.6 进制转换练习

• 使用MATLAB回答以下问题：十进制数123、123.123、321、321.123、223、322、333和333.3的二进制表示是什么？
• 为什么123与123.123、321与321.123、333与333.3的二进制表示相同？

1.7 字符输出练习

编写一个脚本，接受一个输入数字（整数），并在命令窗口中按以下顺序打印出字符：

A
BCD
EFGHI
JKLMNOP
QRSTUVWXY

1.8 矩阵分解练习

• 使用练习11中矩阵[A]和[B]的数值，评估QR、LU、LQ、Cholesky、Schur和奇异值分解。解释为什么[A]和[B]的某些分解（矩阵因式分解）无法计算。
• 使用画廊矩阵函数创建5x5的Hilbert矩阵，并使用Chol_decoposition.m脚本计算Cholesky分解。编辑脚本（Chol_decoposition.m），使其仅计算Cholesky分解的下三角矩阵。
• 使用画廊矩阵函数创建3x3的Riemann矩阵，并计算其QR、LU、LQ、Cholesky、Schur和奇异值分解。
• 执行以下操作：
• 使用画廊矩阵函数创建具有归一化列和指定奇异值的4x4随机矩阵。提示：使用randcolu。
• 计算刚创建矩阵的QR、LU、LQ和分解。
• 执行以下操作：
• 创建一个5x5的随机矩阵，其随机整数元素在1到13的范围内，命名为A_mat。
• 使用Krylov矩阵画廊创建一个5x5的Krylov矩阵，命名为K_mat。
• 使用逻辑运算（A_mat ≥ K_mat）和A_mat与K_mat的逐元素矩阵乘法创建一个名为Logic_A的5x5逻辑值矩阵。
• 创建以下10x10矩阵：

AaA =
    17    24     1     8    15    17    24     1     8    15
    23     5     7    14    16    23     5     7    14    16
     4     6    13    20    22     4     6    13    20    22
    10    12    19    21     3    10    12    19    21     3
    11    18    25     2     9    11    18    25     2     9
    17    24     1     8    15    17    24     1     8    15
    23     5     7    14    16    23     5     7    14    16
     4     6    13    20    22     4     6    13    20    22
    10    12    19    21     3    10    12    19    21     3
    11    18    25     2     9    11    18    25     2     9

提示：使用magic()和repmat()。

2. 常微分方程

2.1 常微分方程定义

许多具有工程应用的建模问题可以用常微分方程（ODEs）来表述。微分方程有几种不同的定义，最简单的一种是“微分方程是包含导数（普通导数或偏导数）的任何方程”。从这个定义可以推导出两种类型的微分方程：常微分方程（ODEs）和偏微分方程（PDEs）。ODEs包含一种类型的导数或一个独立变量，而PDEs则包含两个或更多的导数或独立变量。

例如，一阶ODEs可以表示为：
[\frac{dy}{dx} = f(y, x)]
其中，(y(x))是因变量，其值取决于自变量(x)的值。另一个ODEs的好例子是牛顿第二运动定律，其公式为：
[ma = \frac{dp}{dt} = m\frac{dv}{dt} = F(t, v)]
其中，(F(t, v))是力，是时间(t)和速度(v)的函数；(\frac{dv}{dt})是移动物体的速度变化率（加速度）；(m)是移动物体的质量；(a)是移动物体的加速度；(p)是动量；(\frac{dp}{dt})是其导数。牛顿第二定律的这个公式也可以改写为：
[m\frac{d}{dt}(\frac{dx}{dt}) = m\frac{d^2x}{dt^2} = F(t, x, \frac{dx}{dt})]
其中，移动物体位移(x)的导数(\frac{dx}{dt})是速度(v)。换句话说，速度是移动物体位移(x(t))随时间的变化率。这可以用下面的流程图来可视化：

graph LR
    A[时间t] --> B[位移x]
    B --> C[速度v=dx/dt]
    C --> D[加速度a=dv/dt]
    D --> E[力F(t,v)]
    E --> D

2.2 ODEs分类

ODE相关问题有两种分类：
- 初值问题（IVPs） ：(\ddot{x} = -3\dot{x})，初始条件为(x(0) = 3)，(\dot{x}(0) = 1)。
- 边值问题（BVPs） ：(\ddot{x} = -3\dot{x})，边界条件为(x(0) = 3)，(x(2) = 1.50)。

初值问题（IVPs）由ODEs和未知函数在解域内给定点的指定值（称为初始条件）定义。在ODEs的IVP中，可能有唯一解、无解或多个解。根据定义，ODEs的IVP可以显式或隐式定义。大多数IVP是显式定义的。我们先从显式定义的IVP开始，然后再讨论隐式定义的IVP。此外，根据解的类型（解值在解搜索空间中的变化方式），IVP分为刚性和非刚性问题。而且，ODEs分为线性和非线性两类，并进一步分为齐次和非齐次两组。

ODE类型	示例
刚性ODEs	(\dot{x} = -0.04x + 10^4yz)，(\dot{y} = 0.04x - 10^4yz - 3\times10^7y^2)，(\dot{z} = 3\times10^7y^2)，(t \in [0, 40])
非刚性ODEs	(\dot{y} + 2y = 2t)，(\ddot{w} + 5\dot{w} = 0)
线性ODEs	(\frac{dv}{dt} = -v + 9.81 + 0.198x)，(\ddot{x} + 3\dot{x} + 5x = 0)
非线性ODEs	(\frac{dv}{dt} = -v^2 + 9.81 - 0.198x)，(\ddot{x} + 3\dot{x}^2 + 5x = 0)
齐次ODEs	(\dot{y} + 2y = 0)，(\ddot{x} + 3\dot{x} + 5x = 0)
非齐次ODEs	(\dot{y} + 2y = e^t)，(\ddot{x} + 3\dot{x} + 5x = 2\sin(t))

2.3 ODEs应用示例

• 生物有机体的无约束增长 ：这是一个指数增长问题，描述了生物有机体（如细菌）的无约束增长。这种行为也可以描述房地产或投资价值、流行社交网站的会员增加、零售业务的增长、电气系统的正反馈和化学反应的产生。问题由以下一阶ODE表述：
  [\frac{dy}{dt} = \mu y]
  解为：(y(t) = y_0e^{\mu t})
• 放射性衰变 ：这是指指数衰变，描述了自然界和工程中的许多现象，如放射性衰变、反应堆中化学物质的冲洗、电容器的放电和河流中物质的分解。它用以下一阶ODE表示：
  [\frac{dy}{dt} = -\mu y]
  解为：(y(t) = y_0e^{-\mu t})
• 牛顿第二定律 ：使用牛顿第二定律，下落物体的运动用以下方程表示：
  [m\frac{d^2y}{dt^2} = mg - \gamma\frac{dy}{dt}]
  这是一个二阶ODE，解的形式为：
  [y(t) = C_1e^{-\frac{\gamma}{m}t} + \frac{mg}{\gamma}(1 - e^{-\frac{\gamma}{m}t})]
  其中，(m)是下落物体的质量；(g)是重力加速度；(\gamma)是下落物体的空气阻力系数。三个参数(m)、(g)和(\gamma)是常数，两个参数随时间变化：(\frac{d^2y}{dt^2})（加速度）和(\frac{dy}{dt})（速度）。在下落物体的解中，(C_1)和(C_2)是任意数，取决于初始条件。换句话说，可以根据下落物体的初始条件来计算它们。

2.4 解析方法

符号数学工具箱（或MuPAD笔记本）有几个函数能够评估许多可解析求解的ODEs的解析解。有两个命令（内置函数）可以用来评估某些ODEs的解析解：dsolve和ilaplace/laplace。

2.4.1 DSOLVE

在MATLAB中，计算任何给定ODE的解析（或一般）解的第一个ODE求解工具是dsolve。它可以使用以下一般语法：

Solution = dsolve(equation)
Solution = dsolve(equation, conditions)
Solution = dsolve(equation, conditions, Name, Value)
[y1,...,yN] = dsolve(equations)
[y1,...,yN] = dsolve(equations, conditions)
[y1,...,yN] = dsolve(equations, conditions, Name, Value)

示例1：使用dsolve
给定一个ODE (\dot{y} + 2y^2t = 0)，注意未指定初始和边界条件。

>> y_solution=dsolve('Dy=-2*y^2*t')
Y_solution=-1/(C3-t^2)

注意，(C3)由给定ODE的初始或边界条件定义。
在MATLAB 2012及以后的版本中，也可以使用以下命令语法来解决给定问题：

>> syms y(t); y_sol=dsolve(diff(y) == -2*y^2*t)
 y_sol =
               0
 -1/(- t^2 + C3)

示例2：从dsolve绘图
给定一个ODE (\dot{y} + 2y^2t = 0)，初始条件为(y(0) = 0.5)：

>> Solution=dsolve('Dy=-2*y^2*t', 'y(0)=0.5')
Solution =
 1/(t^2 + 2)

在MATLAB 2012及以后的版本中，替代命令语法如下：

>> syms y(t); Solution=dsolve(diff(y) == -2*y^2*t, y(0)==0.5)
Solution =
 1/(t^2 + 2)

解析解是符号形式的，因此可以使用ezplot进行绘图：

>> ezplot(Solution)

解析解（方程）的数值可以通过向量化（参数化）符号形式（解）来计算：

>> ysol=vectorize(solution)
ysol =
1./(t.^2 + 2)
>> t=(-5:.1:5); ysol_values=eval(ysol);

示例3：添加未指定参数
给定(\dot{y} + ky^2t = 0)，(y(0) = 0.5)，注意这个练习有一个未指定的参数(k)。

>> syms k
>> solution=dsolve('Dy=-k*y^2*t', 'y(0)=0.5')
solution =
1/((k*t^2)/2 + 2)

替代命令语法如下：

>> syms y(t) k; solution=dsolve(diff(y) == -k*y^2*t, y(0)==0.5)
solution =
1/((k*t^2)/2 + 2)

注意：dsolve中的选项需要根据问题类型进行适当设置。在MATLAB 2008 - 2010或更早版本中，应将IgnoreAnalyticConstraints设置为none以获得所有可能的解。例如：

solution=dsolve('Dy=-k*y^2*t', 'y(0)=0.5', ... 'IgnoreAnalyticConstraints', 'none')

对于MATLAB 2012或更高版本，应将IgnoreAnalyticConstraints设置为false以获得所有参数值的所有可能正确答案。否则，由于其预代数简化，dsolve可能会输出不正确的答案。例如：

solution=dsolve(diff(y)==-k*y^2*t, y(0)==0.5, ... 'IgnoreAnalyticConstraints', false)

2.4.2 二阶ODEs和ODE系统

许多过程和现象可以用二阶微分方程来表示。例如，弹簧 - 质量系统的简谐运动、有加速度的物体运动（牛顿第二定律）、阻尼振动、电阻 - 电容 - 电感电路中的电流流动等等。一般来说，二阶ODEs有两种不同的形式：齐次（方程8.4）和非齐次（方程8.5）。
[\ddot{y} + p(x)\dot{y} + q(x)y = 0]
[\ddot{y} + p(x)\dot{y} + q(x)y = g(x)]
注意，方程8.4中的齐次ODEs总是有一个平凡解，即(y(x) = 0)，它满足方程8.4中的给定条件。关于独立函数(p(x))、(q(x))、(g(x))，ODEs可以是线性或非线性的。在某些情况下，独立函数(p(x))、(q(x))、(g(x))可以是常数值或非常数值。

示例1：dsolve求解二阶ODE
给定一个ODE (\frac{d^2u}{dt^2} + 100u = 2.5\sin(10t))，(u(0) = 0)，(\frac{du}{dt}(0) = 0)：

usol=dsolve('D2u+100*u=2.5*sin(10*t)', 'u(0)=0', 'Du(0)=0');
pretty(usol)
%% 替代语法
syms u(t)
Du = diff(u);
u(t) = dsolve(diff(u, 2)==2.5*sin(10*t)-100*u, u(0)==0, Du(0) == 0);
pretty(u(t))

执行这两个短脚本/命令的输出是：

3 sin(10 t)   sin(30 t)
  ----------- - --------- -
      320          320
               / t   sin(20 t) \
     cos(10 t) | - - --------- |
               \ 8      160    /

示例2：ODE系统
给定一个ODE系统：
[\begin{cases}
\dot{y_1} = y_2 \
\dot{y_2} = -y_1 - 0.125y_2
\end{cases}]
初始条件为(y_1(0) = 1)，(y_2(0) = 0)。

%% 用dsolve求解两个一阶ODE系统
yt=dsolve('Dy1=y2', 'Dy2=-y1-0.125*y2','y1(0)=1', 'y2(0)=0');
pretty(yt.y1)
pretty(yt.y2)
%% 替代语法
syms y1(t) y2(t)
z=dsolve(diff(y1,1)==y2, diff(y2,1)==(-y1-0.125*y2), y1(0)==1, y2(0)==0);
pretty(z.y1), pretty(z.y2)

在命令窗口中显示的问题的计算解析解是：

     /    1/2   \             /    1/2   \
     | 255    t |      1/2    | 255    t |
  cos| -------- |   255    sin| -------- |
     \    16    /             \    16    /
  --------------- + ----------------------
          1/16                    1/16
    exp(t)              255 exp(t)
                 /    1/2   \
          1/2    | 255    t |
    16 255    sin| -------- |
                 \    16    /
  - -------------------------
                   1/16
         255 exp(t)

计算得到的解析解为：
[y_1(t) = \frac{16}{16}\cos(\frac{\sqrt{255}}{16}t) + \frac{\sqrt{255}}{16}\sin(\frac{\sqrt{255}}{16}t)]
[y_2(t) = -\frac{16\sqrt{255}}{16}\sin(\frac{\sqrt{255}}{16}t)]

示例3：dsolve无法求解的解
给定一个二阶ODE：(2\ddot{y} + 3\dot{y}^3 - \cos(100t)|y| - 2 = 0)，(y(0) = 0)，(\dot{y}(0) = 0)。

>> Solution_dsolve=dsolve('2*D2y+3*(Dy^3)-cos(100*t)*abs(y)-2=0', ... 
'y(0)=0','Dy(0) =0');
Warning: Explicit solution could not be found.
> In dsolve at 194

这是一种处理无法用dsolve解析求解的ODEs的好方法。

2.5 Laplace变换

具有常系数的线性常微分方程的解可以使用Laplace变换来评估。Laplace变换在求解微分方程中的一个最重要特征是，变换后的方程是一个代数方程，它将用于定义给定微分方程的解。一般来说，Laplace变换应用于求解微分方程可以用以下方式表述。

考虑(y_n(x) = f(t))的(n)阶导数。(y_n(x))的Laplace变换如下：
[\mathcal{L}{\frac{d^n y}{dt^n}} = s^n Y(s) - s^{n - 1}y(0) - s^{n - 2}y’(0) - \cdots - sy^{(n - 2)}(0) - y^{(n - 1)}(0) = F(s)]
或者
[s^n Y(s) - \sum_{i = 1}^{n - 1}s^{n - i - 1}y^{(i)}(0) = F(s)]
在上述方程中，如果我们代入(t = 0)时的初始条件常数值，即(y(0) = a_0)；(y’(0) = a_1)；(y’‘(0) = a_2)；(\cdots)；(y^{n - 2}(0) = a_{n - 2})；(y^{n - 1}(0) = a_{n - 1})。
现在，我们可以将上述表达式重写为：
[\mathcal{L}{\frac{d^n y}{dt^n}} = s^n Y(s) - s^{n - 1}a_0 - s^{n - 2}a_1 - s^{n - 3}a_2 - \cdots - sa_{n - 2} - a_{n - 1} = F(s)]
随后，我们首先求解(Y(s))，然后对(Y(s))取逆Laplace变换，得到(n)阶微分方程的解(y(t))。

应用Laplace和逆Laplace变换来确定具有常系数的微分方程解的一般步骤如下：
1. 对给定方程的两边取Laplace变换。
2. 用(F(s))和其他未知数表示(Y(s))。
3. 对(Y(s))取逆Laplace变换，得到解(y(t))。

2.6 数值方法

由于并非所有的ODEs都能通过解析方法求解，或者解析求解的成本过高，因此数值方法在求解ODEs中具有重要地位。常见的数值方法有欧拉法（向前、向后、改进）、Heun法、中点规则、龙格 - 库塔法、龙格 - 库塔 - 吉尔法、亚当斯 - 巴什福思法、米尔恩法、亚当斯 - 莫尔顿法、泰勒级数法和梯形规则法等。这些方法可分为显式和隐式两类，且根据求解方式又可分为单步方法和多步方法。

方法类型	特点	示例方法
单步方法	仅参考一个先前点及其导数来确定当前值	欧拉法
多步方法	保留并使用先前步骤的值，以获得更高阶的步长	龙格 - 库塔法、亚当斯 - 巴什福思法等

在使用数值方法求解初值问题（IVPs）时，通常从初始点（初始条件）开始，然后在时间上向前迈出一步（等步长或变步长）来计算后续的数值解。所有这些方法都根据其在计算特定类型ODEs的IVP数值解时的准确性和效率（如计算时间和资源消耗）进行评估。

2.7 数值方法的实现

以下以欧拉法为例，介绍如何实现数值方法来求解ODEs。欧拉法是一种简单的单步方法，其基本思想是使用当前点的导数来估计下一个点的值。

考虑一阶ODE (\frac{dy}{dt} = f(t, y))，初始条件为 (y(t_0) = y_0)。欧拉法的迭代公式为：
[y_{n + 1} = y_n + hf(t_n, y_n)]
其中，(h) 是步长，(t_{n + 1} = t_n + h)。

以下是使用欧拉法求解 (\frac{dy}{dt} = -2y^2t)，初始条件 (y(0) = 0.5) 的MATLAB代码示例：

% 定义ODE函数
f = @(t, y) -2 * y^2 * t;

% 初始条件
t0 = 0;
y0 = 0.5;

% 时间范围和步长
t_end = 5;
h = 0.01;

% 初始化时间和y值
t = t0:h:t_end;
y = zeros(size(t));
y(1) = y0;

% 欧拉法迭代
for i = 1:length(t) - 1
    y(i + 1) = y(i) + h * f(t(i), y(i));
end

% 绘制结果
plot(t, y);
xlabel('时间t');
ylabel('y值');
title('欧拉法求解dy/dt = -2y^2t');

上述代码的执行流程如下：
1. 定义ODE函数 (f(t, y))。
2. 设置初始条件 (t_0) 和 (y_0)。
3. 确定时间范围和步长 (h)。
4. 初始化时间向量 (t) 和 (y) 值向量。
5. 使用欧拉法的迭代公式进行迭代计算。
6. 绘制结果。

2.8 不同方法的比较

不同的数值方法在准确性、稳定性和计算效率方面存在差异。一般来说，多步方法通常比单步方法更准确，但计算复杂度也更高。例如，龙格 - 库塔法是一种常用的多步方法，它通过在每个时间步内计算多个斜率来提高准确性。

以下是一个比较欧拉法和四阶龙格 - 库塔法求解同一ODE的示例：

% 定义ODE函数
f = @(t, y) -2 * y^2 * t;

% 初始条件
t0 = 0;
y0 = 0.5;

% 时间范围和步长
t_end = 5;
h = 0.01;

% 初始化时间和y值
t = t0:h:t_end;
y_euler = zeros(size(t));
y_rk4 = zeros(size(t));
y_euler(1) = y0;
y_rk4(1) = y0;

% 欧拉法迭代
for i = 1:length(t) - 1
    y_euler(i + 1) = y_euler(i) + h * f(t(i), y_euler(i));
end

% 四阶龙格 - 库塔法迭代
for i = 1:length(t) - 1
    k1 = h * f(t(i), y_rk4(i));
    k2 = h * f(t(i) + h/2, y_rk4(i) + k1/2);
    k3 = h * f(t(i) + h/2, y_rk4(i) + k2/2);
    k4 = h * f(t(i) + h, y_rk4(i) + k3);
    y_rk4(i + 1) = y_rk4(i) + (k1 + 2*k2 + 2*k3 + k4) / 6;
end

% 绘制结果
figure;
plot(t, y_euler, 'b', 'DisplayName', '欧拉法');
hold on;
plot(t, y_rk4, 'r', 'DisplayName', '四阶龙格 - 库塔法');
xlabel('时间t');
ylabel('y值');
title('欧拉法和四阶龙格 - 库塔法对比');
legend;

从上述代码可以看出，四阶龙格 - 库塔法在每个时间步内计算了四个斜率 (k_1)、(k_2)、(k_3) 和 (k_4)，然后使用加权平均来更新 (y) 值，从而提高了准确性。

2.9 总结

本文介绍了MATLAB中线性代数和常微分方程的相关知识，包括矩阵操作、方程求解、多项式求解、数组创建、声音播放、进制转换、字符输出、矩阵分解等线性代数练习，以及常微分方程的定义、分类、应用示例、解析方法和数值方法。

在常微分方程部分，我们了解了ODEs的基本概念和分类，以及一些常见的解析求解方法（如dsolve和Laplace变换）和数值求解方法（如欧拉法和龙格 - 库塔法）。通过具体的示例和代码，展示了如何使用MATLAB来解决这些问题。

在实际应用中，需要根据具体问题的特点选择合适的方法。对于可解析求解的ODEs，优先考虑使用解析方法，以获得精确的解；对于无法解析求解的ODEs，则需要使用数值方法来获得近似解。同时，不同的数值方法在准确性和效率上存在差异，需要根据实际需求进行选择。

以下是本文内容的流程图总结：

graph LR
    A[线性代数练习] --> B[矩阵操作]
    A --> C[方程求解]
    A --> D[多项式求解]
    A --> E[数组创建]
    A --> F[声音播放]
    A --> G[进制转换]
    A --> H[字符输出]
    A --> I[矩阵分解]
    J[常微分方程] --> K[定义与分类]
    J --> L[应用示例]
    J --> M[解析方法]
    J --> N[数值方法]
    M --> O[dsolve]
    M --> P[Laplace变换]
    N --> Q[欧拉法]
    N --> R[龙格 - 库塔法]

希望本文能够帮助读者更好地理解和应用MATLAB解决线性代数和常微分方程问题。在实际使用中，还可以进一步探索MATLAB的其他功能和工具，以提高解决问题的能力。