MatLab4数值计算一

1、MatLab数学建模22授课：曾景峰（江西环境工程职业学院）第四讲数值计算符号数学工具箱 符号表达式的运算nu meric符号到数值的转换p retty显示悦目的符号输出subs替代子表达式sym建立符号矩阵或表达式symadd符号加法symdiv符号除法symmul符号乘法symop符号运算sympow符号表达式的幕运算symrat有理近似symsub符号减法symvar求符号变量符号表达式的简化collect expand factor simp le sim plify symsum合并同类项展开因式求解最简形式 简化和级数符号多项式char poly hornernu mde npo

2、Iy2sym sym2 poly特征多项式嵌套多项式表示分子或分母的提取多项式向量到符号的转换 符号到多项式向量的转换符号微积分diffint微分积分jorda n taylor约当标准形 泰勒级数展开符号可变精度算术digits设置可变精度vpa可变精度计算求解符号方程compose函数的复合dsolve微分方程的求解fin verse函数逆lin solve齐次线性方程组的求解Solve代数方程的求解符号线性代数char ploy特征多项式determ矩阵行列式的值eige nsys特征值和特征向量in verse矩阵逆jorda n约当标准形lin solve齐次线性方程组的解trans

3、pose矩阵的转置方程求解求解单个代数方程MATLA具有求解符号表达式的工具,如果表达式不是一个方程式（不含等 号），则在求解之前函数solve将表达式置成等于0。 solve( a*xA2+b*x+c )% solve for the roots of the equti onans=1/2/a*(-b+(bA2-4*a*c)A1/2)1/2/a*(-b-(bA2-4*a*c)A1/2)结果是符号向量，其元素是方程的2个解。如果想对非缺省x变量求解，solve 必须指定变量。 solve( a*x2+b*x+c ,b ) % solve for bans=-(a*x2+c)/x带有等号的符号

4、方程也可以求解。 f=solve( cos(x)=sin(x) ) % solve for x f=1/4*pi t=solve( tan(2*x)=sin(x) )t= 0acos(1/2+1/2*3A(1/2)acos(1/2=1/2*3A(1/2)并得到数值解。 numeric (f) ans=0.7854 numeric(t) ans=00 + 0.8314i1.9455注意在求解周期函数方程时，有无穷多的解。在这种情况下， solve 对解的 搜索范围限制在接近于零的有限范围，并返回非唯一的解的子集。如果不能求得符号解，就计算可变精度解。 x=solve( exp(x)=tan(x)

5、 ) x=1.306326940423079代数方程组求解 的未知数求解 n 个方程。可以同时求解若干代数方程，语句 solve （s1 ，s2， ，sn） 对缺省变量求 解n个方程，语句 solve（s1 , s2, . , sn, v1 , v2, . , vn ）对n个v1 , v2， .vnsolve (f)解符号方程式 f 。solve(f1,fn）解由f1,fn组成的联立方程式。我们先定义以下的方程式：eq1 = x-3=4; % 注意也可写成 eq1=x-7eq2 = x2-x-6=0; %注意也可写成eq2=x*2-x-6eq3 = x2+2*x+4=0;eq4 = 3*x+2

6、*y-z=10;eq5 = -x+3*y+2*z=5;eq6 = x-y-z=-1;solve(eq1) ans=solve(eq2)ans=3,-2 %原方程式有二个根3, -2solve(eq3) ans=注意实根和虚根的表示式-1+i*3A(1/2),-1-i*3A(1/2) % solve(eq4,eq5,eq6) % 解三个联立方程式ans= x = -2, y = 5, z = -6如何处理中小学典型的代数问题？黛安娜(Diane)想去看电影，她从小猪存钱罐倒出硬币并清点，她发现:10美分的硬币数加上5美分的硬币总数的一半等于25美分的硬币数。1美分的硬币数比5美分、10美分以及2

7、5美分的硬币总数多10。25美分和10美分的硬币总数等于1美分的硬币数加上1/4的5美分的硬币*硬币数多1。25美分的硬币数和1美分的硬币数比5美分的硬币数加上8倍的10美分的如果电影票价为3.00美元，爆米花为1.00美元，糖棒为50美分，她有足够的 钱去买这三样东西？首先，根据以上给出的信息列出一组线性方程, 美分，5美分，10美分，和25美分的硬币数假如p, n, d和q分别表示1n + pd+=q p = n +d+q10 q+d2q + p = n +8d 1然后，建立MATLAB号方程并对变量求解。 eq1= d+(n+p)/2=q eq2= p=n+d+q-10 eq3= q+d

8、 二p+n/4 eq4= q+p=n+8*d-1 pennies , nickles , dimes, quarters=solve(equ1 p , n, d, q )penni es=16ni ckles=8dimes=3quarters=15,equ2，equ3, equ4,所以，黛安娜有16枚1美分的硬币，8枚 5美分的硬币, 枚25美分的硬币，这就意味着3枚10美分的硬币，15 mon ey=01*16+.05*8+.10*3+.25*15mon ey=4.6100她就有足够的钱去买电影票，爆米花和糖棒并剩余 11美分。【例】求解二元函数方程组Jf1（x, ysin（Xy） =0的零

9、点。lf2（x,y） = cos（x + y） =0沪生x-y平面上网点坐标(0)从三维坐标初步观察两函数图形相交情况 x=-2:0.05:2;y=x;X, Y=meshgrid(x,y);F1=si n(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), shadi ng interp, hold off22y-2-20、图 563-0两函数的三维相交图(1)在某区域观察两函数0等位线的交点情况cle

10、ar;x=-2:0.5:2;y=x;X,Y=meshgrid(x,y);F1=si n(X-Y);F2=cos(X+Y);沪生x-y平面上网点坐标v=-0.2, 0, 0.2;%旨定三个等位值，是为了更可靠地判断0等位线的存在。con tour(X,Y,F1,v)hold on,con tour(X,Y,F2,v),hold off%画 F1的三条等位线。%画 F2的三条等位线。图5.6.3-1 两个二元函数0等位线的交点图(2) 从图形获取零点的初始近似值 在图563-1中，用ginput获取两个函数0等位线(即三线组中间那条线)交点的坐标。%在图上取两个点的坐标x0,y0=g inp ut

11、 (2);dis p( x0,y0)-0.7926 -0.78430.79260.7843(3) 利用fsolve求精确解。以求(0.7926,7843 )附近的解为例。本例直接用字符串表达被解函数。注意：在此，自变量必须写成x(1), x(2)。 假如写成xy(1), xy(2)，指令运行将出错。%fun=si n(x(1)-x (2),cos(x(1)+x(2);xy=fsolve(fu n,x0(2),y0 (2)%xy =0.78540.7854(4) 检验fxy1=si n(xy(1)-xy(2);fxy2=cos(xy(1)+xy (2);dis p(fxy1,fxy2)1.0e-

12、006 *-0.09940.20191说明指令可用以下任何一组指令取代。(A)内联函数形式指令%项x必须有。fun=i niin e(s in (x(1)-x (2), cos(x(1)+x (2), x);xy=fsolve(fu n,x0(2), y0(2);(B) M函数文件形式及指令先用如下fun.m表示被解函数(并在搜索路径上)fun .mfunction ff=fun(x)ff(1)=si n(x(1)-x (2);ff(2)=cos(x(1)+x (2);然后运行指令 xy=fsolve(fun,x0(2),y0(2)。第四步检验中的结果表明：所找零点处的函数值小于10，是一个十

13、分接近零的小数。该精度由 op tio ns.ToIFu n 控制。op tio ns.ToIFu n 的缺省值是 1.0000e-006。它可以用下列指令看到op tio ns=op timset(fsolve);op tio ns.Tol Funans =1.0000e-006线性方程求解a= 7 2 1-215-2-2-21113b=4 7 -1 0x=ab0.49790.14450.0629-0.0813单个微分方程常微分方程有时很难求解，MATLAB供了功能强大的工具，可以帮助求解微 分方程。函数dsovle计算常微分方程的符号解。因为我们要求解微分方程，就需 要用一种方法将微分包含

14、在表达式中。 些不同，用字母D来表示求微分，D2, 程。任何D后所跟的字母为因变量。MATLA解常微分方程式的语法是所以， dsovle 句法与大多数其它函数有一D3等等表示重复求微分，并以此来设定方dsolve(equation,condition)equation代表常微分方程式即y =g(x, y)，且须 以Dy代表一阶微分项y 表二阶微分项 y， condition 则为初始条件。，其中D2y代方程d2y/dx2= 0用符号表达式D2y=(来表示。独立变量可以指定或由symvar规则选定为缺省。例如，一阶方程dy/dx=1+y2的通解为: dsolve( Dy=1+yA2 ) % fi

15、nd the gen eral solutionans=-tan (-X+C1)其中，C1是积分常数。求解初值y(0)=1的同一个方程就可产生:,y(0)=1 ) % add an initial，y(0)=1 ， v ) % find dsolveC Dy=1+yA2 con diti ony=ta n(x+1/4* pi)独立变量可用如下形式指定： dsolve( Dy=1+y2 soluti on tody/dvans=ta n(v+1/4* pi)让我们举一个二阶微分方程的例子，该方程有两个初始条件:2=cos(2x)-y 业(0)=0 y(0)=1 dx2dx y=dsolve( D

16、2y=cos(2*x)-y ， Dy(0)=0 ， y(0)=1 ) y=-2/3*cos(x)A2+1/3+4/3*cos(x) y=si mp le(y) % y looks like it can be simp lified y=-1/3*cos(2*x)+4/3*cos(x)通常，要求解的微分方程含有一阶以上的项，并以下述的形式表示:乌-2 业-3y=0dx dx通解为: y=dsolve( D2y-2Dy-3*y=0 ) y=C1*ex p(-x)+C2*ex p(3*x)加上初始条件：y(0)=0和y(i)=i可得到：,y(o)=o , y(i)=i ) y=solve( D2y

17、-2Dy-3*y=0 y=1/(ex p( -1)-ex p(3)*ex p(-x)-1/(ex p( -1)-ex p( 3)*ex p(3*x) y=si mp le(y)% this looks like a can didateforsim plificatio ny=-(ex p(-x)-ex p(3*x)/(ex p(3)-ex p(-1) P retty(y) % p retty it up exp(-x)-ex p(3 x)exp(3) -exp(-l)现在来绘制感兴趣的区域内的结果。 ezplot(y ，-6 2)例:假设有以下三个一阶常微分方程式和其初始条件y =3x2,

18、y(2)=0.5 y =2x cos(y) 2, y(0)=0.25 CfJ y =3y+ex p( 2x), y(0)=3对应上述常微分方程式的符号运算式为:soln_1 = dsolve(Dy = 3*x2,y(2)=0.5) ans=XA3-7.500000000000000ez plot(soln_1,2,4) %看看这个函数的长相soln_2 = dsolve(Dy = 2*x*cos(y)2,y(0) = pi/4) ans=ata n(xA2+1) soln_3 = dsolve(Dy = 3*y + exp (2*x), y(0) = 3) ans= -ex p(2*x)+4*

19、ex p(3*x)微分方程组函数dsolve也可同时处理若干个微分方程式，下面有两个线性一阶方程。f =3f+4gdx包=-4f+3g dx通解为: f，g=dsolve ( Df=3*f+4*g ，Dg=-4*f+3*g )f=C1*ex p(3*x)*si n(4*x)+C2*ex p(3*x)*cos(4*x)g=-C2*ex p(3*x)*si n(4*x)+C1*ex p(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=ex p(3*x

20、)*si n(4*x)g=ex p(3*x)*cos(4*x)微分和积分微分和积分是微积分学研究和应用的核心，并广泛地用在许多工程学科。MATLA符号工具能帮助解决许多这类问题。微分符号表达式的微分以四种形式利用函数 diff ： f= a*x3+x2-b*x-cexpressionf=a*x3+x2-b*x-c% define a symbolic diff (f) % differentiate with respect to the default variable xans=3*a*xA2+2*x-b diff(f ， a ) % differentiate with respect

21、to aans=x3 diff(f ， 2) % differentiate twice with respect to xans=6*a*x+2 diff(fto aans=0a ，2) % differentiatetwice with respect函数diff也可对数组进行运算。如果F是符号向量或数组，diff(F)对数组内 的各个元素进行微分。 F=sym( a*x symbolic arrayF= a*x c*x3 ，b*x2 ；c*x3d*s ) % create a diff(F) %,b*x2d*sdifferentiate the element with respect

22、to xans= a，2*b*x03*c*x2，注意函数diff也用在MATLAB计算数值向量或矩阵的数值差分。对于一个数 值向量或矩阵M diff(M)计算M(2: m,: )-M(1 : m-1，:)的数值差分，如下 所示： m=(1 : 8)42) % create a vectorM=1 4 9 16 25 36 49 64 diff(M) % find the differences between elements ans=3 5 7 9 11 13 15如果diff的表达式或可变参量是数值，MATLA就非常巧妙地计算其数值差分；如果参量是符号字符串或变量，MATLA对其表达式进行

23、微分。积分积分函数int(f)，其中f是一符号表达式，它力图求出另一符号表达式F使diff(F)=f 。正如从研究微分学所了解的，积分比微分复杂得多。积分或逆求导 不一定是以封闭形式存在， 或许存在但软件也许找不到， 或者软件可明显地求解， 但超过内存或时间限制。当MATLA不能找到逆导数时，它将返回未经计算的命令。 in t( log(x)/ex 9(x2) ) % atte mpt to i ntegrate ans=log(x)/ex p(xA2)同微分一样，积分函数有多种形式。形式 int(f) 相对于缺省的独立变量求逆 导数；形式(f， s )相对于符号变量S积分；形式int(f ，

24、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) % i

25、ntegrate 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)-si n(2*x) int(f ， m ， n ) % integrate with respect to x from m to nans=-1/2*cos(s+2* n)+1/2*cos(s+2*m)正如函数diff 样，积分函数int对符号数组的每一个元素进行运算。 F=sym( a*x， b*x2 ； c*x3，d*s ) % crea

26、te a symbolic arrayF=a*x ，b*x2c*x3 ，d*s diff(F) % ubtegrate the array eleme nts with res pect to x ans=1/2*a*x2，1/3*b*xA31/4*c*x4， d*s*xdiff函数用以演算一函数的微分项，相关的函数语法有下列4个:diff(f)传回f对预设独立变数的一次微分值传回f对独立变数t的一次微分值diff(f, n)传回f对预设独立变数的n次微分值diff(f,t,n)传回 f 对独立变数 t 的 n 次微分值先定义下列三个方程式，接著再演算其微分项：S1 = 6*xA3-4*xA2

27、+b*x-5;S2 = sin(a);S3 = (1 -护3”(1 + 护4);diff(S1) ans= 18*x2-8*x+b diff(S1,2) ans= 36*x-8 diff(S1,b) ans=diff(S2) ans= cos(a) diff(S3) ans= -3*护2心+护4)-4*(1-护3”(1+护4)八2*护3 simplify(diff(S3) ans=护2*(-3+护4-4*圳(1+护4)八2这个函数要找出一符号式 F 使得不存在的int 函数用以演算一函数的积分项， diff(F)=f 。如果积分式的解析式 (analytical form, closed fo

28、rm) 话或是MATLA无法找到，则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*x3-4*x2+b*x-5;S2 = sin(a);S3 = sqrt(x);int(S1)ans=3/2*x4-4/3*x3+

29、1/2*b*x2-5*xint(S2)ans=-cos(a)int(S3)ans=2/3*xA(3/2) i nt(S3,a,b) ans= 2/3*bA(3/2)- 2/3*aA(3/2) i nt(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.074 1数值积分先考虑一个积分式的数学式如下:其中a, b分别为这个积分式的上限及下限，f(x)则为要积分的函树。要求解上述 的积分式，必须设定a, b和f(x)。以MATLAB的积分函数求解的过程亦同， 也要定义f(x)及设定a,b，还须设定在区间a,b之间离散 点(discretized points) 数目，剩下的工作就是选择精度不同的积分法来求解了。梯形法MATLAB提供最简单的积分函数是梯形法trapz,我们先说明梯形法语法 trapz(x,y)，其中x,y分别代