matlab在函数的求解方法 (1)

百度文库- 让每个人平等地提升自我！MatLab & 数学建模
第四讲数值计算
符号表达式的运算
nu meric 符号到数值的转
换
pre tty 显示悦目的符号输出
su
bs
替代子表达式
sy m 建立符号矩阵或表达式
sy
madd
符号加法
sy
mdiv
符号除法
sy
mmul
符号乘法
sy
mop
符号运算
sy mpow 符号表达式的幂运算
sy
mrat
有理近似
sy
msub
符号减法
sy
mvar
求符号变量
符号表达式的
简化
col lect
合并同类项
ex pand 展开
fac
tor
因式
si
mple
求解最
简形式
si
mplify
简化
sy
msum
和级数
符号多项式
char
poly
特征多项式
horn
er
嵌套多项式表示
num
den
分子或分母的提取
poly
2sym
多项式向量到符号
的转换
sym
2poly
符号到多项式向量
的转换
符号微积分
d
iff
微分
i
nt
积分
j
ordan
约当标
准形
t
aylor
泰勒级
数展开
符号可变精度算术
digits 设置可变精度
vpa 可变精度计算
求解符号方程compose 函数的复合dsolve 微分方程的求解finverse 函数逆linsolve 齐次线性方程组的求解solve 代数方程的求解
符号线性代数charploy 特征多项式determ 矩阵行列式的值eigensys 特征值和特征向量inverse 矩阵逆jordan 约当标准形linsolve 齐次线性方程组的解transpose 矩阵的转置
一、方程求解
求解单个代数方程
MATLAB具有求解符号表达式的工具,如果表达式不是一个方程式(不含等号)，则在求解之前函数solve将表达式置成等于0。

>> solve( ' a*x^2+b*x+c ' ) % solve for the roots of the eqution ans=
[1/2/a*(-b+(b^2-4*a*c)^1/2)]
[1/2/a*(-b-(b^2-4*a*c)^1/2)]
结果是符号向量，其元素是方程的2个解。

如果想对非缺省x变量求解，solve 必须指定变量。

>> solve( ' a*x^2+b*x+c ' ， ' b ' ) % solve for b
ans=
-(a*x^2+c)/x
带有等号的符号方程也可以求解。

>> f=solve( ' cos(x)=sin(x) ' ) % solve for x
f=
1/4*pi
>> t=solve( ' tan(2*x)=sin(x) ' )
t=
[ 0]
[acos(1/2+1/2*3^(1/2))]
[acos(1/2=1/2*3^(1/2))]
并得到数值解。

>> numeric(f)
ans=
0.7854
>> numeric(t)
ans=
0 + 0.8314i
1.9455
注意在求解周期函数方程时，有无穷多的解。

在这种情况下，solve对解的搜索范围限制在接近于零的有限范围，并返回非唯一的解的子集。

如果不能求得符号解，就计算可变精度解。

>> x=solve( ' exp(x)=tan(x) ' )
x=
1.306326940423079
代数方程组求解
可以同时求解若干代数方程，语句solve(s1，s2，.....，sn)对缺省变量求解n个方程，语句solve(s1，s2，...，sn，' v1，v2，...，vn ')对n个' v1，v2，...vn '的未知数求解n个方程。

solve(f) 解符号方程式f。

solve(f1,…,fn) 解由f1,…,fn组成的联立方程式。

我们先定义以下的方程式：
>>eq1 = 'x-3=4'; % 注意也可写成'eq1=x-7'
>>eq2 = 'x*2-x-6=0'; % 注意也可写成'eq2=x*2-x-6'
>>eq3 = 'x2+2*x+4=0';
>>eq4 = '3*x+2*y-z=10';
>>eq5 = '-x+3*y+2*z=5';
>>eq6 = 'x-y-z=-1';
>>solve(eq1)
ans=
7
>>solve(eq2)
ans=
[[3],[-2]]' % 原方程式有二个根3, -2
>>solve(eq3)
ans=
[[-1+i*3^(1/2)],[-1-i*3^(1/2)]]' % 注意实根和虚根的表示式
>>solve(eq4,eq5,eq6) % 解三个联立方程式
ans=
x = -2, y = 5, z = -6
如何处理中小学典型的代数问题？
黛安娜(Diane)想去看电影，她从小猪存钱罐倒出硬币并清点，她发现：
•10美分的硬币数加上5美分的硬币总数的一半等于25美分的硬币数。

•1美分的硬币数比5美分、10美分以及25美分的硬币总数多10。

•25美分和10美分的硬币总数等于1美分的硬币数加上1/4的5美分的硬币数
•25美分的硬币数和1美分的硬币数比5美分的硬币数加上8倍的10美分的硬币数多1。

如果电影票价为3.00美元，爆米花为1.00美元，糖棒为50美分，她有足够的钱去买这三样东西？
首先，根据以上给出的信息列出一组线性方程，假如p，n，d和q分别表示1美分，5美分，10美分，和25美分的硬币数
d
n p
q p n d q q d p
n
q p n d
+
+
==++-+=+
+=+-
2
10
4
81
然后，建立MATLAB符号方程并对变量求解。

>> eq1= ' d+(n+p)/2=q ' ；
>> eq2= ' p=n+d+q-10 ' ；
>> eq3= ' q+d=p+n/4 ' ；
>> eq4= ' q+p=n+8*d-1 ' ；
>>[pennies ，nickles ，dimes ，quarters]=solve(equ1，equ2，equ3，equ4，' p ，n ，d ，q ' ) pennies= 16 nickles= 8 dimes= 3
quarters= 15
所以，黛安娜有16枚1美分的硬币，8枚5美分的硬币，3枚10美分的硬币，15枚25美分的硬币，这就意味着
>> money=.01*16+.05*8+.10*3+.25*15 money= 4.6100
她就有足够的钱去买电影票，爆米花和糖棒并剩余11美分。

【例】求解二元函数方程组⎩⎨⎧=+==-=0)cos(),(0
)sin(),(2
1y x y x f y x y x f 的零点。

（0）从三维坐标初步观察两函数图形相交情况
x=-2:0.05:2;y=x;[X,Y]=meshgrid(x,y); %产生x-y 平面上网点坐标 F1=sin(X-Y);F2=cos(X+Y); F0=zeros(size(X)); surf(X,Y,F1),
xlabel('x'),ylabel('y'), view([-31,62]),hold on, surf(X,Y,F2),surf(X,Y,F0), shading interp, %(间隔补齐) hold off
图 5.6.3-0 两函数的三维相交图
（1）在某区域观察两函数0等位线的交点情况 clear;
x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y); %产生x-y 平面上网点坐标 F1=sin(X-Y);F2=cos(X+Y); v=[-0.2, 0, 0.2]; %指定三个等位值，是为了更可靠地判断0等位线的存在。

contour(X,Y,F1,v) %画F1的三条等位线。

hold on,contour(X,Y,F2,v),hold off %画F2的三条等位线。

-2
-1.5-1-0.500.51 1.52
-2-1.5-1-0.500.511.52
图 5.6.3-1 两个二元函数0等位线的交点图
（2）从图形获取零点的初始近似值
在图5.6.3-1中，用ginput 获取两个函数0等位线（即三线组中间那条线）交点的坐标。

[x0,y0]=ginput(2); %在图上取两个点的坐标 disp([x0,y0])
-0.7926 -0.7843 0.7926 0.7843
（3）利用fsolve 求精确解。

以求（0.7926,7843）附近的解为例。

本例直接用字符串表达被解函数。

注意：在此，自变量必须写成x(1), x(2)。

假如写成xy(1), xy(2)，指令运行将出错。

fun='[sin(x(1)-x(2)),cos(x(1)+x(2))]'; %<12> xy=fsolve(fun,[x0(2),y0(2)])
%<13>
xy =
0.7854 0.7854
（4）检验
fxy1=sin(xy(1)-xy(2));fxy2=cos(xy(1)+xy(2));disp([fxy1,fxy2])
1.0e-006 *
-0.0994 0.2019
〖说明〗
●指令<12><13>可用以下任何一组指令取代。

（A）内联函数形式指令
fun=inline('[sin(x(1)-x(2)), cos(x(1)+x(2))]', 'x'); %项'x'必须有。

xy=fsolve(fun,[x0(2), y0(2)]);
（B）M函数文件形式及指令
先用如下fun.m表示被解函数（并在搜索路径上）
[fun.m]
function ff=fun(x)
ff(1)=sin(x(1)-x(2));
ff(2)=cos(x(1)+x(2));
然后运行指令xy=fsolve('fun',[x0(2),y0(2)]) 。

●第四步检验中的结果表明：所找零点处的函数值小于6
10 ，是一个十分接近
零的小数。

该精度由options.TolFun控制。

options.TolFun的缺省值是
1.0000e-006。

它可以用下列指令看到
options=optimset('fsolve');
options.TolFun
ans =
1.0000e-006
线性方程求解
a= [ 7 2 1 -2
9 15 3 -2
-2 -2 11 5
1 3
2 13]
b=[4 7 -1 0]'
x=a\b
x =
0.4979
0.1445
0.0629
-0.0813
单个微分方程
常微分方程有时很难求解，MATLAB提供了功能强大的工具，可以帮助求解微分方程。

函数dsovle计算常微分方程的符号解。

因为我们要求解微分方程，就需要用一种方法将微分包含在表达式中。

所以，dsovle句法与大多数其它函数有一些不同，用字母D来表示求微分，D2，D3等等表示重复求微分，并以此来设定方程。

任何D后所跟的字母为因变量。

MATLAB解常微分方程式的语法是dsolve('equation','condition')，其中equation代表常微分方程式即y'=g(x,y)，且须以Dy代表一阶微分项y'D2y代表二阶微分项y''，condition则为初始条件。

方程d y dx
/=0用符号表达式D2y=0来表示。

独立变量可以指定或由symvar 22
规则选定为缺省。

例如，一阶方程dy/dx=1+y2的通解为：
>> dsolve( ' Dy=1+y^2 ' ) % find the general solution
ans=
-tan(-x+C1)
其中，C1是积分常数。

求解初值y(0)=1的同一个方程就可产生：
>> dsolve(' Dy=1+y^2 '，' y(0)=1 ') % add an initial
condition
y=
tan(x+1/4*pi)
独立变量可用如下形式指定：
>> dsolve(' Dy=1+y^2 '，' y(0)=1 '，' v ') % find
solution to dy/dv
ans=
tan(v+1/4*pi)
让我们举一个二阶微分方程的例子，该方程有两个初始条件：
d y dx
22=cos(2x)-y dy dx (0)=0 y(0)=1
>> y=dsolve(' D2y=cos(2*x)-y '，' Dy(0)=0 '，' y(0)=1 ') y=
-2/3*cos(x)^2+1/3+4/3*cos(x)
>> y=simple(y) % y looks like it can be simplified y=
-1/3*cos(2*x)+4/3*cos(x)
通常，要求解的微分方程含有一阶以上的项，并以下述的形式表示：
d y dx
22
-2dy
dx -3y=0 通解为：
>> y=solve( 'D2y-2Dy-3*y=0 ')
y=
C1*exp(-x)+C2*exp(3*x)
加上初始条件：y(0)=0和y(1)=1可得到：
>> y=solve( ' D2y-2Dy-3*y=0 ' ， ' y(0)=0，y(1)=1 ' ) y=
1/(exp(-1)-exp(3))*exp(-x)-1/(exp(-1)-exp(3))*exp(3*x)
>> y=simple(y) % this looks like a candidate for simplification y=
-(exp(-x)-exp(3*x))/(exp(3)-exp(-1))
>> pretty(y) % pretty it up exp(-x)-exp(3 x) - --------------------- exp(3) -exp(-1)
现在来绘制感兴趣的区域内的结果。

>> ezplot(y，[-6 2])
-(exp(-x)-exp(3*x))/(exp(3)-exp(-1))
5
-5
-10
-15
-6-5-4-3-2-1012
x
例:假设有以下三个一阶常微分方程式和其初始条件
y'=3x2, y(2)=0.5
y'=2.x.cos(y)2, y (0)=0.25
y'=3y+exp(2x), y(0)=3
对应上述常微分方程式的符号运算式为：
>>soln_1 = dsolve('Dy = 3*x^2','y(2)=0.5')
ans=
x^3-7.500000000000000
>>ezplot(soln_1,[2,4]) % 看看这个函数的长相
>>soln_2 = dsolve('Dy = 2*x*cos(y)^2','y(0) = pi/4')
ans=
atan(x^2+1)
>>soln_3 = dsolve('Dy = 3*y + exp(2*x)',' y(0) = 3')
ans=
-exp(2*x)+4*exp(3*x)
微分方程组
函数dsolve 也可同时处理若干个微分方程式，下面有两个线性一阶方程。

dx
df =3f+4g dg dx =-4f+3g
通解为：
>> [f ，g]=dsolve ( ' Df=3*f+4*g ' ， ' Dg=-4*f+3*g ' )
f=
C1*exp(3*x)*sin(4*x)+C2*exp(3*x)*cos(4*x)
g=
-C2*exp(3*x)*sin(4*x)+C1*exp(3*x)*cos(4*x)
加上初始条件：f(0)=0和g(0)=1，我们可以得到：
>> [f ，g]=dsolve( ' Df=3*f+4*g ' ， ' Dg=-4*f+3*g ' ， ' f(0)=0，
g(0)=1 ' )
f=
exp(3*x)*sin(4*x)
g=
exp(3*x)*cos(4*x)
微分和积分
微分和积分是微积分学研究和应用的核心，并广泛地用在许多工程学科。

MATLAB符号工具能帮助解决许多这类问题。

微分
符号表达式的微分以四种形式利用函数diff：
>> f= ' a*x^3+x^2-b*x-c ' % define a symbolic expression
f=
a*x^3+x^2-b*x-c
>> diff(f) % differentiate with respect to the default variable x
ans=
3*a*x^2+2*x-b
>> diff(f，'a ') % differentiate with respect to a
ans=
x^3
>> diff(f，2) % differentiate twice with respect to x
ans=
6*a*x+2
>> diff(f，' a '，2) % differentiate twice with respect to a
ans=
函数diff也可对数组进行运算。

如果F是符号向量或数组，diff(F)对数组内的各个元素进行微分。

>> F=sym(' [a*x， b*x^2； c*x^3， d*s] ') % create a
symbolic array
F=
[ a*x， b*x^2]
[c*x^3， d*s]
>> diff(F) % differentiate the element with respect to x
ans=
[ a，2*b*x]
[3*c*x^2， 0]
注意函数diff也用在MATLAB，计算数值向量或矩阵的数值差分。

对于一个数值向量或矩阵M，diff(M)计算M(2： m，： )-M(1： m-1，： )的数值差分，如下所示：
>> m=[(1： 8).^2)] % create a vector
M=
1 4 9 16 25 36 49 64
>> diff(M) % find the differences between elements
ans=
3 5 7 9 11 13 15
如果diff的表达式或可变参量是数值，MATLAB就非常巧妙地计算其数值差分；如果参量是符号字符串或变量，MATLAB就对其表达式进行微分。

积分
积分函数int(f)，其中f是一符号表达式，它力图求出另一符号表达式F使diff(F)=f。

正如从研究微分学所了解的，积分比微分复杂得多。

积分或逆求导不一定是以封闭形式存在，或许存在但软件也许找不到，或者软件可明显地求解，但超过内存或时间限制。

当MATLAB不能找到逆导数时，它将返回未经计算的命令。

>> int( ' log(x)/exp(x^2) ' ) % attempt to integrate
ans=
log(x)/exp(x^2)
同微分一样，积分函数有多种形式。

形式int(f)相对于缺省的独立变量求逆导数；形式(f，' s ')相对于符号变量s积分；形式int(f，a，b)和int(f，' s '，a，b)，a，b是数值，求解符号表达式从a到b的定积分；形式int(f，' m ' ，' n ')和形式int(f，' s '，' m '，' n ')，其中m，n是符号变量，求解符号表达式从m到n的定积分。

>> f=' sin(s+2*x) ' % crate a symbolic function
f=
sin(s+2*x)
>> int(f) % integrate with respect to x
ans=
-1/2*cos(s+2*x)
>> int(f，' s ') % integrate with respect to s
ans=
-cos(s+2*x)
>> int(f，pi/2，pi) % integrate with respect to x from
π/2 toπ
ans=
-cos(x)
>> int(f，' s '，pi/2，pi) % integrate with respect to
s from π/2 to π
ans=
cos(2*x)-sin(2*x)
>> int(f，' m '，' n ') % integrate with respect to x
from m to n
ans=
-1/2*cos(s+2*n)+1/2*cos(s+2*m)
正如函数diff一样，积分函数int对符号数组的每一个元素进行运算。

>> F=sym( ' [a*x，b*x^2；c*x^3，d*s] ' ) % create a symbolic
array
F=
[ a*x，b*x^2]
[c*x^3， d*s]
>> diff(F) % ubtegrate the array elements with respect to x
ans=
[1/2*a*x^2，1/3*b*x^3]
[1/4*c*x^4， d*s*x]
diff函数用以演算一函数的微分项，相关的函数语法有下列4个：
diff(f) 传回f对预设独立变数的一次微分值
diff(f,'t')传回f对独立变数t的一次微分值
diff(f,n)传回f对预设独立变数的n次微分值
diff(f,'t',n)传回f对独立变数t的n次微分值
先定义下列三个方程式，接著再演算其微分项：
>>S1 = '6*x^3-4*x^2+b*x-5';
>>S2 = 'sin(a)';
>>S3 = '(1 - t^3)/(1 + t^4)';
>>diff(S1)
ans=
18*x^2-8*x+b
>>diff(S1,2)
ans=
36*x-8
>>diff(S1,'b')
ans=
x
>>diff(S2)
ans=
cos(a)
>>diff(S3)
ans=
-3*t^2/(1+t^4)-4*(1-t^3)/(1+t^4)^2*t^3
>>simplify(diff(S3))
ans=
t^2*(-3+t^4-4*t)/(1+t^4)^2
int函数用以演算一函数的积分项，这个函数要找出一符号式 F 使得
diff(F)=f。

如果积分式的解析式 (analytical form, closed form) 不存在的
话或是MATLAB无法找到，则 int 传回原输入的符号式。

相关的函数语法有下列4个：
int(f)传回f对预设独立变数的积分值
int(f,'t')传回f对独立变数t的积分值
int(f,a,b)传回f对预设独立变数的积分值，积分区间为[a,b]，a和b为数值式
int(f,'t',a,b)传回f对独立变数t的积分值，积分区间为[a,b]，a和b为数值式
int(f,'m','n')传回f对预设变数的积分值，积分区间为[m,n]，m和n为符号式
我们示范几个例子：
>>S1 = '6*x^3-4*x^2+b*x-5';
>>S2 = 'sin(a)';
>>S3 = 'sqrt(x)';
>>int(S1)
ans=
3/2*x^4-4/3*x^3+1/2*b*x^2-5*x
>>int(S2)
ans=
-cos(a)
>>int(S3)
ans=
2/3*x^(3/2)
>>int(S3,'a','b')
ans=
2/3*b^(3/2)- 2/3*a^(3/2)
>>int(S3,0.5,0.6)
ans=
2/25*15^(1/2)-1/6*2^(1/2)
>>numeric(int(S3,0.5,0.6)) % 使用numeric函数可以计算积分的数值
ans=
0.0741
数值积分
先考虑一个积分式的数学式如下：
其中a, b分别为这个积分式的上限及下限，f(x) 则为要积分的函树。

要求解上述的积分式，必须设定a, b和f(x)。

以MATLAB 的积分函数求解的过程亦同，也要定义f(x) 及设定a,b，还须设定在区间[a,b] 之间离散点(discretized points) 数目，剩下的工作就是选择精度不同的积分法来求解了。

梯形法
MATLAB提供最简单的积分函数是梯形法trapz，我们先说明梯形法语法
trapz(x,y)，其中x,y分别代表数目相同的阵列或矩阵，而y与x的关系可以由是一函数型态（如y=sin(x)）或是不以函数描述的离散型态。

我们看一简单积分式
以下为MATLAB 的程式
>> x=0:pi/100:pi;
>> y=sin(x);
>> k=trapz(x,y)
k =
1.9998
二次函数法
MATLAB 另外提供二种积分函数，它们分别是辛普森法quad和牛顿-康兹法quad8。

三种方法的精确度由低而高，分别为trapz, quad, quad8。

由于这二种方法依据的积分法不同于梯形法，因此它们的语法就和trapz不同；其语法为quad('function',a,b)（quad8语法相同），其中function是一已定义函数的名称（如sin, cos, sqrt, log等），而a, b是积分的下限和上限。

和trapz
比较，quad, quad8不同之处在于这二者类似解析式的积分式，只须设定上下限及定义要积分的函数；而trapz则是针对离散点型态的数据做积分。

我们看一简单积分式
以下为MATLAB 的程式
>> a=0; b=0.5;
>> kq=quad('sqrt',a,b)
kq =
0.2357
>> kq8=quad8('sqrt',a,b)
kq8 =
0.2357
再来看一个较复杂的积分式
>> x=-1:0.17:2;
>> y=humps(x);
>> area=trapz(x,y)
area =
25.9174
>> x=-1:0.07:2;
>> y=humps(x);
>> area=trapz(x,y)
area =
26.6243
>> area=quad('hump',-1,2) area =
26.3450
>> area=quad8('hump',-1,2) area =
26.3450
符号表达式画图
在许多的场合，将表达式可视化是有利的。

MATLAB 提供了函数ezplot 来完成该任务。

>> y=' 16*x^2+64*x+96 ' % expression to plot
y=
16*x^2+64*x+96
>> ezplot(y)
图 符号函数16*x^2+64*x+96 (-2π≤x ≤2π)
ezplot 绘制了定义域为-2π≤x ≤2π的给定符号函数，并相应地调整了y 轴比例，还加了网格栅和标志。

在这个例子中，我们感兴趣的时间是从0到6。

让我们再试一下，并指定时间范围，
>> ezplot(y ，[0 6]) % plot y for 0≤x ≤6 。

-6-4-20
246-1000
-800
-600
-400
-200
200
x -16*x^2+64*x+96
-16*x^2+64*x+96
150
100
50
-50
-100
0123456
x
图符号函数16*x^2+64*x+96 0≤x≤6
21。