MATLAB 之级数计算和方程符号求解

这篇具有很好参考价值的文章主要介绍了MATLAB 之级数计算和方程符号求解。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

一、级数

级数是表示函数，研究函数性质以及进行数值计算的一种工具，特别是可以利用收敛的无穷级数来逼近一些无理数，使它们的求值变得更方便。

1. 级数符号求和

前面曾讨论过有限级数求和的函数 sum，sum 处理的级数是以一个向量形式表示的，并且只能是有穷级数，对于无穷级数求和，sum 是无能为力的。
求无穷级数的和需要符号表达式求和函数 symsum，其调用格式如下：

	symsum(s,v,n,m)

其中， $s$ 表示一个级数的通项，是一个符号表达式。 $v$ 是求和变量， $v$ 省略时使用系统的默认变量。 $n$ 和 $m$ 是求和的开始项和末项。
例如，我们求下列定积分。
（1） $s_{1}=1+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+\cdots +\frac{1}{n^{2}}+\cdots$
（2） $s_{2}=1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+\cdots +(-1)^{n+1}\frac{1}{n}+\cdots$
（3） $s_{3}=x+2x^{2}+3x^{3}+\cdots +nx^{n}+\cdots$
（4） $s_{4}=1+4+9+16+\cdots +10000$
程序如下：

>> n=sym('n');
>> x=sym('x');
>> s1=symsum(1/n^2,n,1,inf)            	%求s1
 
s1 =
 
pi^2/6
 
>> s2=symsum((-1)^(n+1)/n,1,inf)    	%求s2。未指定求和变量，默认为n
 
s2 =
 
log(2)
 
>> s3=symsum (n*x^n,n,1,inf)           	%求s3。此处的求和变量n不能省略   
 
s3 =
 
piecewise(abs(x) < 1, x/(x - 1)^2)
 
>> s4=symsum(n^2,1,100)              	%求s4。计算有限级数的和
 
s4 =
 
338350

其中，s3 的结果表示 $\left | x \right | < 1$ 时，级数和 s3 $=\frac{x}{(x-1)^2}$ 。

2. 函数的泰勒级数

泰勒（Taylor）级数将一个任意函数表示为一个幂级数，并且，在许多情况下，只需要取幂级数的前有限项来表示该函数，这对于大多数工程应用问题来说，精度已经足够。
MATLAB 提供了 taylor 函数将函数展开为幂级数，其调用格式如下：

	taylor (f,v,a)
	taylor (f, v,a, Name, Value)

该函数将函数 $f$ 按变量 $v$ 展开为泰勒级数 $v$ 的默认值与 symvar 函数指示的默认变量相同。参数 $a$ 指定将函数 $f$ 在自变量 $v = a$ 处展开， $a$ 的默认值是 0。
第二种格式用于运算时设置相关选项，Name 和 Value 成对使用，Name 为选项，Value 为 Name 的值。Name 有以下 3 个可取字符串。
（1） 'ExpansionPoint'：指定展开点，对应值为标量或向量。未设置时，展开点为 0。
（2） 'Order'：指定截断阶，对应值为一个正整数。未设置时，截断阶为 6，即展开式的最高阶为 5。
（3） 'OrderMode'：指定展开式采用绝对阶或相对阶，对应值为 ‘Absolute’（绝对阶）或 ‘Relative’（相对阶）。未设置时取 ‘Absolute’。
例如，我们求函数在指定点的泰勒级数展开式。
（1）求 $\sqrt{1-2x+x^{3}}-\sqrt[3]{1-3x+x^{2}}$ 的 5 阶泰勒级数展开式。
（2）将 $\sqrt{1-2x+x^{3}}-\sqrt[3]{1-3x+x^{2}}$ 在 $x = 1$ 处按 5 次多项式展开。
程序如下：

>> x=sym('x');
>> f1=sqrt(1-2*x+x^3)-(1-3*x+x^2)^(1/3);
>> f2=(1+x+x^2)/(1-x+x^2);
>> taylor(f1)					%求（1）
 
ans =
 
(239*x^5)/72 + (119*x^4)/72 + x^3 + x^2/6
 
>> taylor(f2,x,1,'Order',6)		%求（2）。展开到x-1的5次幂时应选择x=6
 
ans =
 
2*(x - 1)^3 - 2*(x - 1)^2 - 2*(x - 1)^5 + 3

二、方程符号求解

除了常规的代数方程以及常微分方程数值求解的方法，MATLAB 还提供了 solve 和 dsolve 函数，可以用符号运算求解代数方程和常微分方程。

1. 代数方程符号求解

代数方程是指未涉及微积分运算的方程，相对比较简单。在 MATLAB 中，求解符号表达式表示的代数方程可由函数 solve 实现，其调用格式如下：
① solve(s)：求解符号表达式 $s$ 的代数方程，求解变量为默认变量。
② solve(s,v)：求解符号表达式 $s$ 的代数方程，求解变量为 $v$ 。
③ solve(s1,s2,...,sn],[v1,v2,...,vn])：求解符号表达式 $s_{1},s_{2},\dots ,s_{n}$ 组成的代数方程组，求解变量分别为 $v_{1},v_{2},\dots ,v_{n}$ 。
例如，我们求解下列方程。
（1） $\frac{1}{x+2}+\frac{4x}{x^{2}-4}=1+\frac{2}{x-2}$
（2） $x-\sqrt[3]{x^{3}-4x-7}=1$
（3） $2\sin (3x-\frac{\pi}{4})=1$
（4） $x+xe^{x}-10=0$
程序如下：

>> syms x
>> f=1/(x+2)+4*x/(x^2-4)==1+2/(x-2)				%定义（1）中函数
 
f =
 
1/(x + 2) + (4*x)/(x^2 - 4) == 2/(x - 2) + 1
 
>> x=solve(f)									%解方程（1）
 
x =
 
1
 
>> syms x
>> f=x-(x^3-4*x-7)^(1/3)==1						%定义（2）中方程
 
f =
 
x - (x^3 - 4*x - 7)^(1/3) == 1

>> x=solve(f)									%解方程（2）
 
x =
 
3
 
>> syms x
>> f=2*sin(3*x-pi/4)==1							%定义（3）中方程

f =
 
2*sin(3*x - pi/4) == 1
 
>> x=solve(f)									%解方程（3）
 
x =
 
 (5*pi)/36
(13*pi)/36
 
>> syms x
>> f=x+x*exp(x)-10								%定义（4）中方程
 
f =
 
x + x*exp(x) - 10
 
>> x=solve(f,x)									%解方程（4）。仅标出方程的左端
 
x =
 
1.6335061701558463841931651789789

例如，我们求解下列方程。
（1） $\left\{\begin{matrix}\frac{1}{x^{3}}+\frac{1}{y^{3}}=28 \\\frac{1}{x}+\frac{1}{x}=4 \end{matrix}\right.$
（2） $\left\{\begin{matrix}x+y=98 \\\sqrt[3]{x}+\sqrt[3]{y}=2 \end{matrix}\right.$
（3） $\left\{\begin{matrix}u^{3}+v^{3}=98 \\u+v=2 \end{matrix}\right.$
（4） $\left\{\begin{matrix}x^{2}+y^{2}=5 \\2x^{2}-3xy-2y^{2}=0 \end{matrix}\right.$
程序如下：

>> syms x y
>> f1=1/x^3+1/y^3==28	%定义（1）中方程
 
f1 =
 
1/x^3 + 1/y^3 == 28
 
>> f2=1/x+1/y==4	%定义（1）中方程
 
f2 =
 
1/x + 1/y == 4
 
>> [x,y]=solve([f1,f2],[x,y])	%解方程组（1）
 
x =
 
  1
1/3
 
 
y =
 
1/3
  1
 
>> syms x y
>> f1=x+y==98	%定义（2）中方程
 
f1 =
 
x + y == 98
 
>> f2=x^(1/3)+y^(1/3)==2	%定义（2）中方程
 
f2 =
 
x^(1/3) + y^(1/3) == 2
 
>> [x,y]=solve([f1,f2],[x,y])	%解方程组（2）
 
x =
 
Empty sym: 0-by-1
 
 
y =
 
Empty sym: 0-by-1

对方程组（2），MATLAB 给出了无解的结论，显然错误，请看完全与其同构的方程组（3）。输入程序如下：

>> syms u v
>> f1=u^3+v^3==98	%定义（3）中方程
 
f1 =
 
u^3 + v^3 == 98
 
>> f2=u+v==2	%定义（3）中方程
 
f2 =
 
u + v == 2
 
>> [u,v]=solve([f1,f2],[u,v])	%解方程组（3）
 
u =
 
 5
-3
 
 
v =
 
-3
 5
 
>> syms x y
>> f1=x^2+y^2==5	%定义（4）中方程
 
f1 =
 
x^2 + y^2 == 5
 
>> f2=2*x^2-3*x*y-2*y^2==0	%定义（4）中方程
 
f2 =
 
2*x^2 - 3*x*y - 2*y^2 == 0
 
>> [x,y]=solve([f1,f2],[x,y])	%解方程组（4）
 
x =
 
 1
-2
 2
-1
 
 
y =
 
-2
-1
 1
 2

在应用 MATLAB 求解方程组时，应充分发挥其功能。当 MATLAB 给出无解或未找到所期望的解时，不要认为原方程组就真正无解了，还需要用其他方法试探着求解。
如果知道方程组在某点附近有解，不妨用方程组的数值求解函数 fsolve 试探求解，一般能找到所期望的解。

2. 常微分方程符号求解

在 MATLAB 中，用大写字母 D 表示导数。例如，Dy 表示 $y^{'}$ ，D2y 表示 $y "$ ，Dy(0)=5 表示 $y^{'} (0) = 5$ 。D3y+D2y+Dy -x+5=0 表示微分方程 $y " + y " + y^{'} x + 5 = 0$ 。
符号常微分方程求解可以通过函数 dsolve 来实现，其调用格式如下：

	dsolve(e,c,v)

该函数求解常微分方程 $e$ 在初值条件 $c$ 下的特解。参数 $v$ 描述方程中的自变量，省略时按默认原则处理，若没有给出初值条件 $c$ ，则求方程的通解。
dsolve 函数在求常微分方程组时的调用格式如下：

	dsolve(e1,e2,…,en,c1,…,cn,v)

该函数求解常微分方程组 e1,…,en 在初值条件 c1,…,cn 下的特解，若不给出初值条件，则求方程组的通解。 $v$ 给出求解变量，如果没有指定自变量，则采用默认自变量 $t$ 。
例如，我们求下列定积分。
（1）求 $\frac{\mathrm{d} y}{\mathrm{d} x} =\frac{x^{2}+y^{2}}{2x^{2}}$ 的通解。
（2）求 $x^{2}\frac{\mathrm{d} y}{\mathrm{d} x} +2xy=e^{x}$ 的通解。
（3）求 $\frac{\mathrm{d} y}{\mathrm{d} x} =\frac{x^{2}}{1+y^{2}}$ 的特解， $y (2) = 1$ 。
（4）求 $\left\{\begin{matrix}\frac{\mathrm{d} y}{\mathrm{d} x} =4x-2y \\\frac{\mathrm{d} y}{\mathrm{d} x} =2x-y \end{matrix}\right.$ 的通解。
程序如下：

>> y=dsolve('Dy-(x^2+y^2)/x^2/2','x')		%解（1）。方程的右端为0时可以不写
 
y =
 
-x*(1/(C1 + log(x)/2) - 1)
                         x
 
>> y=dsolve('Dy*x^2+2*x*y-exp(x)','x')		%解（2）

y =
 
-(C1 - exp(x))/x^2
 
>> y=dsolve('Dy-x^2/(1+y^2)','y(2)=1','x')	%解（3）
 
y =
 
(((x^3/2 - 2)^2 + 1)^(1/2) + x^3/2 - 2)^(1/3) - 1/(((x^3/2 - 2)^2 + 1)^(1/2) + x^3/2 - 2)^(1/3)
 
>> [x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y','t')	%解方程组（4）
  
x =
 
C1/2 + 2*C2*exp(3*t)
  
y =
 
C1 + C2*exp(3*t)