MATLAB解方程

MATLAB解方程篇分为解方程基础篇，常微分方程(ODE)篇，偏微分方程(PDE)篇。此篇为解方程基础篇，在本文中，我们首先介绍解方程的基础知识，然后将使用MATLAB解初高中所学习的多元多次方程及方程组。

符号计算与数值计算的概念

首先需要理清的概念是符号计算与数值计算，符号计算是凭借一系列恒等式、数学定理，通过推理和演绎，获得解析结果的运算手段。数值计算则是使用计算机求数学问题近似解的方法与过程。

那么这两者有什么区别？

>> 2*pi	% 数值计算
ans = 6.2832
>> sym(2*pi)	% 符号计算
ans = 2*pi

通过实例我们可以体会到二者差异，数值计算是取了近似解(数值解)，通过对 π \pi π(pi)取近似值的方法将结果计算出来；而符号计算则将pi当做一个变量，得到了一个绝对准确但有些许滑稽的结果 2*pi。

那么符号计算就是这样的，建立在数值完全准确表达和推演严格解析的基础之上，所得结果完全准确的计算方法。

而数值计算则是采用某种计算手段，如有限元方法、数值逼近、插值法等而得到的解。通过数值计算得到的解并不完全准确，但数值计算胜在运算速度快而且可以对函数没有解析解的情况进行处理。

在进行符号计算的学习前应当先了解什么是数值计算以及如何进行数值计算：数值计算与MATLAB微积分

符号计算

符号计算又称计算机代数，它可以使用计算机推导数学公式，如对表达式进行因式分解、化简、微分、积分、解代数方程、求解常微分方程等。

为了进行上面所说的因式分解等符号计算，我们必须要首先定义符号。

构建符号变量

符号常量是不含变量的符号表达式，用sym函数来创建；符号变量使用sym和syms函数来创建。

sym构造符号变量和对象。S=sym(s, 参数)执行时，S为所创建的符号对象，参数表示转换后的格式，可以省略，当被转换的s是数值时，参数可以是’d’, ‘f’, ‘e’或’r’ 4种格式，当被转换的’s’是字符串时，参数可以是’integer’, ‘rational’, ‘real’, ‘positive’, 'clear’等格式。

syms s1 s2 ... sn 等同于 s1=sym(‘s1’), s2=sym(‘s2’)…

在计算物理中，我们更多的使用符号变量。例如使用符号变量构建方程(符号表达式)。使用符号变量，我们可以做到求解代数方程。

clc; clear all
syms a b c x y
f = a*x^2+b*x+c;	% 构建方程
solve(f)		% solve为求解方程的函数，其中默认自变量为x
% output: ans = [-(b + (b^2 - 4*a*c)^(1/2))/(2*a), -(b - (b^2 - 4*a*c)^(1/2))/(2*a)]
solve(f,a)	% 求解方程时可以自定义自变量
% output: ans = -(c + b*x)/x^2

符号表达式的变换

MATLAB的符号数学工具箱提供了collect、expand、horner、factor、simplify和simple等函数实现符号表达式的化简。

计算反函数

函数f(x)存在一个反函数，g(f(x))=x，则g和f互为反函数，MATLAB中finverse函数可以用来求符号函数的反函数

g=finverse(f,v)	% 对f(v)按指定自变量v求反函数，当v省略，则对默认的自由符号变量求反函数

% 求符号函数 f=5*sin(x)+y 的反函数
syms x y
f=5*sin(x)+y
g=finverse(f) 	% 对默认自由变量求反函数
% output: g = -asin(1/5*y-1/5*x)
g1=finverse(f,'y')
% output: g1 = -5*sin(x)+y

计算复合函数

MATLAB提供了compose函数可以求出f(x)和g(y)的复合函数f(g(y))

compose(f,g,x,y,z)	% 计算f和g的复合函数，xyz都省略时计算出复合函数为f(g(y))
% xy省略则计算出f(g(z))，都不省略则以x为自由变量计算出f(g(z))并用z代替y

多项式符号表达式

多项式函数如果是包含分子和分母，则可以使用numden函数来进行通分计算，并可以提取分子和分母多项式符号表达式。

[N,D]=numden(s)	% 提取多项式符号表达式s的分子和分母，N和D是通分过后的分子和分母多项式符号表达式

f1=sym(1/(x-1)+1/(x+1)+3);
[N1,D1]=numden(f1)
% output: N1 = 2*x+3*x^2-2
% output: D1 = x^2-2

符号微积分、极限和级数

微积分–微分使用diff函数，积分使用int函数实现

diff(f,t,n)		% 计算f对符号变量t的n阶微分

% 计算符号表达式f=sin(ax)+y^2*cos(x)的微分
syms a x y
f=sin(a*x)+y^2*cos(x);
dfdx=diff(f)	%对默认自由变量x求一阶微分
% output: dfdx = cos(a*x)*a-y^2*sin(x)
dfdy=diff(f,y)	%对符号变量
% output: dfdy = 2*y*cos(x)
dfdy2=diff(f,y,2)	%对符号变量y求二阶微分
% output: dfdy2 = 2*cos(x)

% 计算符号矩阵的一阶微分
f=sym([sqrt(1+x^2) t^x;exp(x) x+y]);
dfdx=diff(f)
% 输出如下：
%		dfdx = 
%		[ x/(x^2 + 1)^(1/2), t^x*log(t)]
%		[            exp(x),          1]

% diff函数还可以用来计算数组各行的差值
a=[1 2 3;1 0 1];
diff(a)		% output: ans = 0 -2 -2

积分–例： f = ∫ 0 2 π ∫ 0 a r 2 sin ⁡ 2 φ d r d φ f=\int_0^{2\pi}\int_0^ar^2\sin^2{\varphi}drd\varphi f=∫02π​∫0a​r2sin2φdrdφ

int(f,t,a,b)		% 计算符号变量t的积分

% 根据微分表达式f计算原函数
syms a r phi
g=r^2*(sin(phi))^2;
f=int(int(g,r,0,a),phi,0,2*pi)		% output: f = (pi*a^3)/3

极限

limt(f)表达式 lim ⁡ x → 0 f ( x ) \lim_{x\rightarrow0}f(x) limx→0​f(x)，limt(f,a)表达式 lim ⁡ x → a f ( x ) \lim_{x\rightarrow a}f(x) limx→a​f(x)，limt(f,t,a)表达式 lim ⁡ t → a f ( x ) \lim_{t\rightarrow a}f(x) limt→a​f(x)，limt(f,x,a,left)表达式 lim ⁡ x → a − f ( x ) \lim_{x\rightarrow a^-}f(x) limx→a−​f(x)。

syms t
f1=exp(-t)*sin(t);
ess=limit(f1,t,inf)	%计算趋向于无穷大的规模
% output: ess = 0

f2=1/t;
limitf2 = limit(f2)		% output: limitf2=NaN

limitf2_l = limit(f2,'t','0','left')		% output: limitf2_l = -Inf

limitf2_r = limit(f2,'t','0','right')		% output: limitf2_r = Inf

级数求和

symsum(s,x,a,b)	% 计算表达式s当x从a到b的级数和

syms x k n
f1=x^k;
limitf1=symsum(f1,k,0,inf)	% 计算无穷级数	output: limitf1 = -1/(x-1)
f2=1/k;
sumf2=symsum(f2,k,1,10)	% 计算前10项级数	output: sumf2 = 7381/2520

Taylor级数 (泰勒级数)

如果函数 f ( x ) f(x) f(x)在点 x 0 x_0 x0​的某一邻域内具有从一阶到 n + 1 n+1 n+1阶的导数，则在该邻域内函数 f ( x ) f(x) f(x)在点 x = x 0 x=x_0 x=x0​时，趋向无穷的幂级数为 f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + … f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\dots f(x)=f(x0​)+f′(x0​)(x−x0​)+2!f′′(x0​)​(x−x0​)2+… 这个级数称为泰勒级数，MATLAB中使用Taylor函数来计算。

taylor(f,x,x0,'Order',n)	%求泰勒级数以符号变量x在x0点展开n项

% 使用taylor函数对符号表达式cos(x)和e^{-t}*sin(t)进行Taylor级数展开
syms x
f1=cos(x);
taylorf1=taylor(f1,x,1,'Order',3)	%计算x=1级数展开前3项
output: taylorf1 = cos(1) - sin(1)*(x-1) - (cos(1)*(x-1)^2)/2

f2=exp(-x)*sin(x); 
taylorf2=taylor(f2)	% output: taylorf2 = -x^5/30 + x^3/3 - x^2 + x	级数展开到x^5

解方程基本案例 解一般方程

依次对一元多次方程 x 2 − x − 2 = 0 x^2-x-2=0 x2−x−2=0、参数方程 a x 2 + b x + c = 0 ax^2+bx+c=0 ax2+bx+c=0、 常微分方程 d d t x = − a x \frac{d}{dt}x=-ax dtd​x=−ax 求解。

clc;clear
syms a b c x y	% 字符型数据--定义
% 解一元多次方程
f=x^2-x-2;
solve(f) 	% 或solve(f,x)，只有一个变量时变量可以不表示
solve(x^2-x == 2) 		% 也可将solve内的式子写为等于0的等式

% 解参数方程
f=a*x^2+b*x+c;
solve(f,x)

% 常微分方程
dsolve('Dx = -a*x')		% output: ans=C1*exp(-a*t)

解方程组

% 解方程组
f=x^2-y;
g=y-x-2;
[a,b]=solve(f,g)

% 解参数方程组
f=a*x^2+y^2;g=x-y-1;
[sx,sy]=solve(f,g,x,y)

符号对象的运算

符号工具箱有三种不同类型的算术运算方式，包括数值型、有理数型和VPA型。

1)数值型 数值型运算是MATLAB的浮点运算，运算速度快但结果不准确，例如2/3就是一般的浮点运算。

2)有理数型 有理数型的运算是精确符号运算，计算时间和占用内存都最大，其结果是非常精确的，例如sym(2/3)是有理数的算术运算。

3)VPA型 VPA型的算术运算是任意精度运算

%% 符号对象的转换
digits(n)		% 设定n位有效位数的精度，默认为32位
S=vpa(s,n)		% 将s按n为有效数字计算得出符号对象S
n=double(s)		% 将符号对象s转换为双精度数值对象n
% 转为为数值型用double，转换为VPA型用vpa，数值型转换为有理数型用sym

digits	% output: Digits = 32, 显示默认精度
q=sym(sqrt(2))		% output: q = 2^(1/2)
q=vpa(q)	% output: 1.4142135623730950488016887242097, 按默认精度计算并显示

digits(15)
p=vpa(pi)		% output: 3.14159265358979
x=vpa(q*p,10)		% output: 4.442882938
n1=double(x)		% output: 4.4429

更多MATLAB知识学习见：计算物理入门