MatLab & 数学建模第四讲数值计算一、方程求解求解单个代数方程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 bans=-(a*x^2+c)/x带有等号的符号方程也可以求解。

>> f=solve( ' cos(x)=sin(x) ' ) % solve for xf=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.8314i1.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美分的硬币数dn pq p n d q q d pnq p n d++==++-+=++=+-210481然后，建立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= 3quarters= 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(),(21y 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的三条等位线。

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

[x0,y0]=ginput(2); %在图上取两个点的坐标disp([x0,y0])-0.7926 -0.78430.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)]) 。

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

该精度由options.TolFun控制。

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

它可以用下列指令看到options=optimset('fsolve');options.TolFunans =1.0000e-006线性方程求解a= [ 7 2 1 -29 15 3 -2-2 -2 11 51 32 13]b=[4 7 -1 0]'x=a\bx =0.49790.14450.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 solutionans=-tan(-x+C1)其中，C1是积分常数。

求解初值y(0)=1的同一个方程就可产生：>> dsolve(' Dy=1+y^2 '，' y(0)=1 ') % add an initialconditiony=tan(x+1/4*pi)独立变量可用如下形式指定：>> dsolve(' Dy=1+y^2 '，' y(0)=1 '，' v ') % findsolution to dy/dvans=tan(v+1/4*pi)让我们举一个二阶微分方程的例子，该方程有两个初始条件：d ydx22=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 dx22-2dydx -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)现在来绘制感兴趣的区域内的结果。