MatLab在函数的求解方法.doc

1 MatLab % 注意也可写成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) % 注意实根和虚根的表示式 4 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 np qpndqqdp n qpnd 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 5 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 1 yxyxf yxyxf （0 0）从三维坐标初步观察两函数图形相交情况）从三维坐标初步观察两函数图形相交情况 x=-2:0.05:2;y=x;X,Y=meshgrid(x,y);x=-2:0.05:2;y=x;X,Y=meshgrid(x,y);% %产生产生 x-yx-y 平面上网点坐标平面上网点坐标 F1=sin(X-Y);F2=cos(X+Y);F1=sin(X-Y);F2=cos(X+Y); F0=zeros(size(X);F0=zeros(size(X); surf(X,Y,F1),surf(X,Y,F1), xlabel(x),ylabel(y),xlabel(x),ylabel(y), view(-31,62),holdview(-31,62),hold on,on, surf(X,Y,F2),surf(X,Y,F0),surf(X,Y,F2),surf(X,Y,F0), shadingshading interp,interp, %(%(间隔补齐间隔补齐) ) holdhold offoff 图 5.6.3-0 两函数的三维相交图 6 （1 1）在某区域观察两函数）在某区域观察两函数 0 0 等位线的交点情况等位线的交点情况 clear;clear; x=-2:0.5:2;y=x;X,Y=meshgrid(x,y);x=-2:0.5:2;y=x;X,Y=meshgrid(x,y);% %产生产生 x-yx-y 平面上网点坐标平面上网点坐标 F1=sin(X-Y);F2=cos(X+Y);F1=sin(X-Y);F2=cos(X+Y); v=-0.2,v=-0.2, 0,0, 0.2;0.2; % %指定三个等位值，是为了更可靠地判断指定三个等位值，是为了更可靠地判断 0 0 等位线的存在。等位线的存在。 contour(X,Y,F1,v)contour(X,Y,F1,v)% %画画 F1F1 的三条等位线。的三条等位线。 holdhold on,contour(X,Y,F2,v),holdon,contour(X,Y,F2,v),hold offoff % %画画 F2F2 的三条等位线。的三条等位线。 05 15 2 0 5 1 5 2 图 5.6.3-1 两个二元函数 0 等位线的交点图 （2 2）从图形获取零点的初始近似值）从图形获取零点的初始近似值 在图 5.6.3-1 中，用 ginput 获取两个函数 0 等位线（即三线组中间那条线） 交点的坐标。 x0,y0=ginput(2);x0,y0=ginput(2); % %在图上取两个点的坐标在图上取两个点的坐标 disp(x0,y0)disp(x0,y0) -0.7926 -0.7843 0.7926 0.7843 （3 3）利用）利用 fsolvefsolve 求精确解。以求（求精确解。以求（0.7926,78430.7926,7843）附近的解为例。）附近的解为例。 本例直接用字符串表达被解函数。注意：在此，自变量必须写成 x(1), x(2)。 假如写成 xy(1), xy(2)，指令运行将出错。 fun=sin(x(1)-x(2),cos(x(1)+x(2);fun=sin(x(1)-x(2),cos(x(1)+x(2);% xy=fsolve(fun,x0(2),y0(2)xy=fsolve(fun,x0(2),y0(2) % xy = 0.7854 0.7854 （4 4）检验）检验 fxy1=sin(xy(1)-xy(2);fxy2=cos(xy(1)+xy(2);disp(fxy1,fxy2)fxy1=sin(xy(1)-xy(2);fxy2=cos(xy(1)+xy(2);disp(fxy1,fxy2) 1.0e-006 * -0.0994 0.2019 说明 指令可用以下任何一组指令取代。 （A）内联函数形式指令 fun=inline(sin(x(1)-x(2),fun=inline(sin(x(1)-x(2), cos(x(1)+x(2),cos(x(1)+x(2), x);x); % %项项xx必须有。必须有。 7 xy=fsolve(fun,x0(2),xy=fsolve(fun,x0(2), y0(2);y0(2); （B）M 函数文件形式及指令 先用如下 fun.m 表示被解函数（并在搜索路径上） fun.mfun.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)xy=fsolve(fun,x0(2),y0(2) 。 第四步检验中的结果表明：所找零点处的函数值小于，是一个十分接 6 10 近零的小数。该精度由 options.TolFunoptions.TolFun 控制。 options.TolFunoptions.TolFun 的缺省值 是 1.0000e-006。它可以用下列指令看到 options=optimset(fsolve);options=optimset(fsolve); options.TolFunoptions.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=ab x = 0.4979 0.1445 0.0629 -0.0813 单个微分方程单个微分方程 8 常微分方程有时很难求解，MATLAB提供了功能强大的工具，可以帮助求解 微分方程。函数dsovledsovle计算常微分方程的符号解。因为我们要求解微分方程， 就需要用一种方法将微分包含在表达式中。所以，dsovledsovle句法与大多数其它函 数有一些不同，用字母D D来表示求微分，D2D2，D3D3等等表示重复求微分，并以此来 设定方程。任何D D后所跟的字母为因变量。 MATLAB解常微分方程式的语法是dsolve(equation,condition)，其中 equation代表常微分方程式即y=g(x,y)，且须 以Dy代表一阶微分项y D2y代 表二阶微分项y ，condition则为初始条件。 方程=0用符号表达式D2y=0来表示。独立变量可以指定或由symvarsymvard ydx 22 / 规则选定为缺省。例如，一阶方程dy/dx=1+y2的通解为： dsolve( Dy=1+y2 ) % find the general solution ans= -tan(-x+C1) 其中，C1C1是积分常数。求解初值y(0)=1的同一个方程就可产生： dsolve( Dy=1+y2 ， y(0)=1 ) % add an initial condition y= tan(x+1/4*pi) 独立变量可用如下形式指定： dsolve( Dy=1+y2 ， y(0)=1 ， v ) % find solution to dy/dv ans= tan(v+1/4*pi) 让我们举一个二阶微分方程的例子，该方程有两个初始条件： d y dx 2 2 =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= 9 -1/3*cos(2*x)+4/3*cos(x) 通常，要求解的微分方程含有一阶以上的项，并以下述的形式表示： d y dx 2 2 -2 dy 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) -6-5-4-3-2-1012 -15 -10 -5 0 5 x -(exp(-x)-exp(3*x)/(exp(3)-exp(-1) 10 例例:假设有以下三个一阶常微分方程式和其初始条件 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*x2,y(2)=0.5) ans= x3-7.500000000000000 ezplot(soln_1,2,4) % 看看这个函数的长相 soln_2 = dsolve(Dy = 2*x*cos(y)2,y(0) = pi/4) ans= atan(x2+1) soln_3 = dsolve(Dy = 3*y + exp(2*x), y(0) = 3) ans= -exp(2*x)+4*exp(3*x) 微分方程组微分方程组 函数dsolvedsolve也可同时处理若干个微分方程式，下面有两个线性一阶方程。 =3f+4g dg dx =-4f+3g dx df 通解为： 11 ff，g=dsolveg=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符号工具能帮助解决许多这类问题。 微分微分 符号表达式的微分以四种形式利用函数diffdiff： f= a*x3+x2-b*x-c % define a symbolic expression f= a*x3+x2-b*x-c diffdiff(f) % differentiate with respect to the default variable x ans= 3*a*x2+2*x-b diff(f，a ) % differentiate with respect to a 12 ans= x3 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= 0 函数diffdiff也可对数组进行运算。如果F是符号向量或数组，diff(F)diff(F)对数组 内的各个元素进行微分。 F=sym( a*x， b*x2； c*x3， d*s ) % create a symbolic array F= a*x， b*x2 c*x3， d*s diff(F) % differentiate the element with respect to x ans= a，2*b*x 3*c*x2， 0 注意函数diffdiff也用在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就对其表达式进行微分。 积分积分 13 积分积分函数int(f)int(f)，其中f是一符号表达式，它力图求出另一符号表达式F使 diff(F)=fdiff(F)=f。正如从研究微分学所了解的，积分比微分复杂得多。积分或逆求导 不一定是以封闭形式存在，或许存在但软件也许找不到，或者软件可明显地求 解，但超过内存或时间限制。当MATLAB不能找到逆导数时，它将返回未经计算 的命令。 int( log(x)/exp(x2) ) % attempt to integrate ans= log(x)/exp(x2) 同微分一样，积分函数有多种形式。形式int(f)int(f)相对于缺省的独立变量求 逆导数；形式(f(f， s s )相对于符号变量s积分；形式int(fint(f，a a，b)b)和int(fint(f， s s ，a a，b)b)，a a，b b是数值，求解符号表达式从a a到b b的定积分；形式int(fint(f， m m ， n n )和形式int(fint(f， s s ， m m ， n n )，其中m m，n n是符号变量，求解 符号表达式从m m到n 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) 14 正如函数diffdiff一样，积分函数intint对符号数组的每一个元素进行运算。 F=sym( a*x，b*x2；c*x3，d*s ) % create a symbolic array F= a*x，b*x2 c*x3， d*s diff(F) % ubtegrate the array elements with respect to x ans= 1/2*a*x2，1/3*b*x3 1/4*c*x4， 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*x3-4*x2+b*x-5; S2 = sin(a); S3 = (1 - t3)/(1 + t4); diff(S1) ans= 18*x2-8*x+b diff(S1,2) ans= 15 36*x-8 diff(S1,b) ans= x diff(S2) ans= cos(a) diff(S3) ans= -3*t2/(1+t4)-4*(1-t3)/(1+t4)2*t3 simplify(diff(S3) ans= t2*(-3+t4-4*t)/(1+t4)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 为符 号式 16 我们示范几个例子： S1 = 6*x3-4*x2+b*x-5; S2 = sin(a); S3 = sqrt(x); int(S1) ans= 3/2*x4-4/3*x3+1/2*b*x2-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 17 数值积分数值积分 先考虑一个积分式的数学式如下： 其中 a, b 分别为这个积分式的上限及下限，f(x) 则为要积分的函树。要求解上 述的积分式，必须设定 a, b 和 f(x)。以 MATLAB 的积分函数求解的过程亦同， 也要定义 f(x) 及设定 a,b，还须设定在区间 a,b 之间离散 点(discretized points) 数目，剩下的工作就是选择精度不同的积分法来求解了。 梯形法梯形法 MATLAB 提供最简单的积分函数是梯形法 trapz，我们先说明梯形法语法 trapz(x,y)，其中 x,yx,y 分别代表数目相同分别代表数目相同 的阵列或矩阵的阵列或矩阵，而 y 与 x 的关系可以由 是一函数型态（如 y=sin(x)）或是不以函数描述的离散型态。 我们看一简单积分式 以下为 MATLAB 的程式 x=0:pi/100:pi; y=sin(x); k=trapz(x,y) k = 1.9998 二次函数法二次函数法 MATLAB 另外提供二种积分函数，它们分别