MATLAB符号数学：微分方程、表达式转换与应用

# MATLAB 符号数学：微分方程、表达式转换与应用 ## 1. 微分方程基础 微分方程包含因变量及其相对于自变量的导数。例如，$\frac{dy}{dt} = y$ 就是一个微分方程，其中 $y$ 是因变量，$t$ 是自变量。在 MATLAB 中，常微分方程的默认自变量是 $t$。 以 $y = e^t$ 为例，其关于 $t$ 的导数为 $\frac{dy}{dt} = e^t$，由于 $y = e^t$，所以该方程也可表示为 $\frac{dy}{dt} = y$。微分方程的解是用自变量 $t$ 表示因变量 $y$ 的表达式，且通常有多个解。如 $y = C_1e^t$ 这一族函数都可由 $\frac{dy}{dt} = y$ 表示。若指定初始条件，如 $y(0) = 1$，则可确定常数 $C_1 = 1$。 再看 $y = t^2$，其导数为 $\frac{dy}{dt} = 2t$，可改写为 $\frac{dy}{dt} = \frac{2y}{t}$，这同样是一个微分方程。 ## 2. 用 dsolve 函数求解微分方程 ### 2.1 一阶微分方程 MATLAB 的符号工具箱中的 `dsolve` 函数可用于求解微分方程。求解一阶微分方程时，首先要定义符号变量。以 $\frac{dy}{dt} = \frac{2y}{t}$ 为例，操作步骤如下： ```matlab % 定义符号变量 syms y(t) % 定义微分方程 eqn = diff(y(t),t)==y(t)/t; % 或者使用另一种语法 eqn = diff(y,t)==y/t; % 求解微分方程 dsolve(eqn) ``` 输出结果为 `C1*t^2`。若指定边界条件，如 $y(-1) = 1$，可在 `dsolve` 函数的第二个输入字段中输入该条件： ```matlab dsolve(eqn,y(-1)==1) ``` 输出结果为 `t^2`，此时 $C_1 = 1$。 ### 2.2 高阶微分方程 高阶微分方程的求解方法与一阶类似。以 $100\frac{d^2x}{dt^2} + 4\frac{dx}{dt} + 20x = 0$ 为例，边界条件为 $x(0) = 1$ 和 $\frac{dx}{dt}(0) = 1$，求解步骤如下： ```matlab % 定义符号变量 syms x(t) Dx = diff(x); % 定义微分方程 eqn = 100*diff(x,t,2)+4*diff(x,t)+20*x==0; % 求解微分方程 dsolve(eqn,x(0)==1,Dx(0)==1) ``` 输出结果为一个复杂的表达式。 ### 2.3 微分方程组 对于微分方程组，如 $\frac{dx}{dt} = y$ 和 $\frac{dy}{dt} = x$，求解步骤如下： ```matlab % 定义符号函数 syms x(t) y(t) % 求解微分方程组 a = dsolve(diff(x)==y,diff(y)==x); % 查看结果 a.y a.x ``` 结果以结构体数组形式返回。可按常规方式添加边界条件，线性微分方程组中每个变量需一个边界条件，若包含高阶导数则需添加更多边界条件。 需要注意的是，并非所有微分方程都能通过 MATLAB 进行符号求解，对于无法解析求解的方程，需采用数值技术。 ## 3. 将符号表达式转换为匿名函数 在使用传统 MATLAB 函数前，先对数学表达式进行符号计算通常很有用。`matlabFunction` 函数可将符号表达式转换为匿名函数。 ### 3.1 简单示例 ```matlab % 定义符号变量 syms x y = cos(x); dy = diff(y); % 将符号变量转换为匿名函数 f = matlabFunction(dy); % 计算 f(2) f(2) ``` 输出结果为 `-0.9093`。 ### 3.2 复杂示例 ```matlab % 定义符号变量 syms x y(x) = (exp(-x)-1)/x; dy = diff(y); % 将符号变量转换为匿名函数 g = matlabFunction(dy); ``` 得到匿名函数 `g = @(x)-exp(-x)./x-1.0./x.^2.*(exp(-x)-1.0)`。 ## 4. MATLAB 符号变量与表达式操作 ### 4.1 创建符号变量与表达式 使用 `syms` 命令可创建符号变量，如： ```matlab syms x syms a b c ``` 可使用这些符号变量创建更复杂的表达式，如二次方程： ```matlab EX = a*x^2 + b*x + c; ``` ### 4.2 符号表达式与方程的区别 方程是等式，表达式不是。例如： ```matlab EQ = n == m/MW; % 定义符号方程 ``` ### 4.3 符号表达式与方程的操作 MATLAB 提供了许多用于操作符号表达式和方程的函数，如下表所示： | 函数名 | 功能 | | ---- | ---- | | `collect` | 合并同类项 | | `diff` | 求符号表达式的导数 | | `double` | 将符号变量转换为双精度浮点数 | | `dsolve` | 求解微分方程 | | `expand` | 展开表达式或方程 | | `ezpolar` | 创建极坐标图 | | `fcontour` | 创建等高线图 | | `fimplicit` | 根据隐式方程创建二维图 | | `fimplicit3` | 根据隐式方程创建三维图 | | `fmesh` | 根据符号表达式创建网格图 | | `fplot` | 绘制符号表达式的图像 | | `fplot3` | 创建三维线图 | | `fsurf` | 根据符号表达式创建表面图 | | `factor` | 分解表达式或方程 | | `int` | 求符号表达式的积分 | | `matlabFunction` | 将符号表达式转换为匿名函数 | | `numden` | 提取表达式的分子和分母 | | `poly2sym` | 将多项式数组转换为符号多项式 | | `separateUnits` | 将变量分离为单位和值的组件 | | `simplify` | 根据内置规则简化表达式或方程 | | `solve` | 求解符号表达式或方程 | | `subs` | 替换符号表达式或方程中的变量 | | `sym` | 创建符号变量、表达式或方程 | | `sym2poly` | 将符号多项式转换为数组 | | `symfun` | 创建符号函数 | | `syms` | 创建符号变量 | | `symunit` | 向工作区添加物理测量单位 | | `symvar` | 确定函数中的符号变量并按字母顺序返回 | | `unitConvert` | 转换单位 | ## 5. 符号代数问题示例 ### 5.1 创建符号表达式 ```matlab % 创建符号变量 syms a b c d x se1 = x^3 - 3*x^2 + x; se2 = sin(x) + tan(x); se3 = (2*x^2 - 3*x - 2)/(x^2 - 5*x); se4 = (x^2 -9)/(x+3); ``` ### 5.2 表达式的运算 ```matlab % 除法运算 se1/se2; % 乘法运算 se3*se4; % 除法运算 se1/x; % 加法运算 se1 + se3; ``` ### 5.3 创建符号函数 ```matlab % 创建符号变量 syms a b c d x fse1(x) = x^3 - 3*x^2 + x; fse2(x) = sin(x) + tan(x); fse3(x) = (2*x^2 - 3*x - 2)/(x^2 - 5*x); fse4(x) = (x^2 -9)/(x+3); ``` ### 5.4 函数的运算 ```matlab % 除法运算 fse1/fse2; % 乘法运算 fse3*fse4; % 除法运算 fse1/x; % 加法运算 fse1 + fse3; ``` ### 5.5 创建符号方程 ```matlab sq1 = x^2 + y^2 == 4; sq2 = 5*x^5 - 4*x^4 + 3*x^3 + 2*x^2 - x == 24; sq3 = sin(a) + cos(b)-x*c == d; sq4 = (x^3 -3*x)/(x-3) == 14; ``` ### 5.6 使用 numden 函数 尝试使用 `numden` 函数提取 `se4`、`fse4` 和 `sq4` 的分子和分母： ```matlab [n1, d1] = numden(se4); [n2, d2] = numden(fse4); [n3, d3] = numden(sq4); ``` 该函数对表达式有效，但对等式无效。 ### 5.7 使用 expand、factor、collect 和 simplify 函数 ```matlab % 对 se1 到 se4 进行操作 expand(se1); factor(se1); collect(se1); simplify(se1); % 对 fse1 到 fse4 进行操作 expand(fse1); factor(fse1); collect(fse1); simplify(fse1); % 对 sq1 到 sq4 进行操作 expand(sq1); factor(sq1); collect(sq1); simplify(sq1); ``` ## 6. 符号求解与 subs 命令的使用 ### 6.1 求解表达式 ```matlab % 求解 se1 到 se4 中的 x solve(se1, x); solve(se2, x); solve(se3, x); solve(se4, x); ``` ### 6.2 求解函数 ```matlab % 求解 fse1 到 fse4 中的 x solve(fse1, x); solve(fse2, x); solve(fse3, x); solve(fse4, x); ``` ### 6.3 求解方程 ```matlab % 求解 sq1 到 sq4 中的 x solve(sq1, x); solve(sq2, x); solve(sq3, x); solve(sq4, x); ``` ### 6.4 求解特定方程中的变量 ```matlab % 求解 sq3 中的 a solve(sq3, a); ``` ### 6.5 摆锤问题 已知摆锤的频率公式为 $f = \frac{1}{2\pi}\sqrt{\frac{mgL}{I}}$，求解长度 $L$： ```matlab % 定义符号变量 syms f m g L I eqn = f == (1/(2*pi))*sqrt((m*g*L)/I); % 求解 L solve(eqn, L); ``` 若已知 $m = 10$ kg，$f = 0.2$ s$^{-1}$，$I = 60$ kg·m$^2$，$g = 9.8$ m/s$^2$，使用 `subs` 函数求解 $L$： ```matlab % 定义符号变量 syms f m g L I eqn = f == (1/(2*pi))*sqrt((m*g*L)/I); % 求解 L L_sol = solve(eqn, L); % 代入数值 subs(L_sol, [m, f, I, g], [10, 0.2, 60, 9.8]); ``` ### 6.6 动能问题 动能公式为 $KE = \frac{1}{2}mV^2$，求解速度 $V$： ```matlab % 定义符号变量 syms KE m V eqn = KE == (1/2)*m*V^2; % 求解 V solve(eqn, V); ``` 计算质量为 2000 lbm、速度为 60 mph 的汽车的动能，并将单位转换为 Btu： ```matlab % 定义符号变量 syms m V KE eqn = KE == (1/2)*m*V^2; % 代入数值 KE_value = subs(eqn, [m, V], [2000, 60]); % 单位转换 % 1 lbf = 32.174 lbm·ft/s^2 % 1 h = 3600 s % 1 mile = 5280 ft % 1 Btu = 778.169 ft·lbf KE_value_Btu = subs(KE_val ```