《代数方程和常微分方程求解》由会员分享，可在线阅读，更多相关《代数方程和常微分方程求解（40页珍藏版）》请在金锄头文库上搜索。
1、第8章 代数方程和常微分方程求解n代数方程是未知数和常数进行有限次代数运算所 组成的方程，它包括有理方程和无理方程。代数 方程 的解称为 的根或零点， 其求解一般是通过代数几何来进行。n微分方程是含有一个或是多个导数的方程。只有 一个自变量及其导数的微分方程称为常微分方程 ；包含两个以上自变量及其偏导数的微分方程称 为偏微分方程。n工程上许多物理规律，设计过程的模拟和评价， 凡是涉及质量和能量运动设计分析的问题，都最 终归结到微分方程。8.1 代数方程求解n8.1.1 代数方程图解法n符号绘图函数fplot()和ezplot()也可以用于图解 法求代数方程的根，它适用于求解维数较少的一 维方程
2、或二维方程组。n对于一维方程图解，其解就是函数曲线与x轴交点 所对应的变量数值。如果有多个交点，则表示该 方程有多个解；如果没有交点，则表示该方程没 有解。n例如，在例5-3使用符号绘图函数绘制代数方程的 图形（图5-3左图）中可见，函数在区间-5,5内 与x轴有3个交点，因此该代数方程该区间内有3个 实根。n对于二维方程组图解，其解就是两条函数曲线的 交点所对应的坐标数值。如果只有1个交点（或切 点），则表示该方程组有1个解；如果有2个交点 ，则表示该方程组有2个解；如果没有交点，则表 示该方程没有解。n例8-1 用图解法求解二维联立方程。na=-2;b=2; % 定义横轴区间nezplot
3、(x2+y2-1.69,a,b); naxis equal; % 控制坐标轴比例相等nhold on;grid on;nezplot(2.4*x3-y+1.5,a,b); nline(a,b,0,0);line(0,0,b,a);nxlabel(bf x);ylabel(bf y);ntitle(bf 二维代数方程组的图解法) ngtext(bf f_1=x2+y2-1.32);ngtext(bf f_2=2.4x3-y+1.5);n8.1.2 代数方程的解析解n求非线性方程或方程组解析解的函数调用格式：n X=solve(fun,x)n其中，fun是符号方程的函数表达式，x是自变量 ，X是解
4、析解。n应当指出，函数solve(fun,x)也可以用于求线性 方组的解析解。n例8-2 求非线性解方程组解析解n n% 二维非线性方程组的解析解nsyms a b x y;nf1=x2-x*y-a;nf2=y2-x*y+b;ndisp( 二维非线性方程组的解析解：)nX,Y=solve(f1,f2,x,y)nM文件运行结果：n二维非线性方程组的解析解：nx =n a/(a-b)(1/2)n -a/(a-b)(1/2)nY =n 1/(a-b)(1/2)*bn -1/(a-b)(1/2)*bn8.1.3 线性方程组的数值解n最简便方法是使用矩阵左除或是矩阵求逆的方法 ，求解线性方程组AX=b。
5、n X= Abn X=inv(A)*bn其中，A是方程组的系数矩阵，b是常数向量，X是 解析解。n例8-3 求线性方程组的数值解n% 线性方程组的数值解nAA=1,1,1,1;1,2,-1,4;2,-3,-1,-5;3,1,2,11;nbb=5;-2;-2;0; % 线性方程组常数向量ndisp( 采用矩阵左除求出线性方程组的解：)nxx=AAbbndisp( 采用矩阵求逆求出线性方程组的解：)nzx=inv(AA)*bbndisp( 计算残量：)nr=AA*zx-bbndisp( 计算残量的模：)nR=norm(r)nM文件运行结果：n采用矩阵左除或矩阵求逆求出线性方程组的解：nxx (zx
6、)= 1.0000n 2.0000n 3.0000n -1.0000n计算残量：nr = 1.0e-014 *n 0.0888n 0.2220n -0.4441n 0.1776n计算残量的模：nR = 5.3475e-015n8.1.4 非线性方程的数值解n1、一维非线性方程n对于一维非线性方程求解，可以看作是单变量的 极小化问题，通过不断缩小搜索区间来逼近一维 问题的真解。因此，可以使用一维非线性方程优 化解函数来求解。其调用格式是：n x,fx,flag=fzero(fun,x0)n其中，输入参数中：fun是非线性方程的函数表达 式；x0是根的初值；n输出参数中：x是非线性方程的数值解；f
7、x是数值 解的函数值；返回参数flag0时，表示求解成功 ，否则求解出现问题。n函数fzero所使用的算法为二分法、secant法和逆 二次插值法的组合。n例8-4 求解一维非线性方程n% 求解单变量x非线性方程nx0=0.1; % 解的初值nxz,fz,flag=fzero(atan(x)+exp(x),x0);ndisp( 求解成功性判断参数：), flagndisp( 非线性方程的解：),xzndisp( 非线性方程解的函数值：),fznM文件运行结果：n求解成功性判断参数：nflag = 1n非线性方程的解：nxz = -0.6066n非线性方程解的函数值：nfz = -1.1102e
8、-016n2、多维非线性方程组n对于多维非线性方程组n使用多维非线性方程组优化解函数求解的格式：n x,fval,flag=fsolve(fun,x0)n其中，输入参数中fun是非线性方程组的向量函数 表达式；x0是根的初值；输出参数中x是非线性方 程（组）的数值解；fval是数值解的函数值；返 回参数flag0时，表示求解成功.n函数fsolve的作用是从根的初值x0开始，以逐渐 减少误差的算法，搜索出满足多维非线性方程组 fun的实根x和对应的函数值fval。n例8-5 求解二维非线性方程组nx0=2;3; % 解的初值n% 定义非线性方程组表达式f和向量xnfun=inline(2*x(
9、1)-x(2)2-exp(-x(1);- x(1)3+x(1)*x(2)-exp(-x(2),x)nxn,fval,flag=fsolve(fun,x0);ndisp( 求解成功性判断参数：), flagnif flag0n disp( 方程组的解成功)nelseif flag0，解成功。8.2 常微分方程求解n求解微分方程必须事先对自变量的某些值规定出 函数或是导数的值。n若在自变量为零的点上，给出初始条件，称为初 值问题，最普遍的自变量是“时间”。例如，弹性 系统的自由振动，若以时间为零来限定位移和速 度，这是一个初值问题。n若在自变量为非零的点上，给出边界条件，称为 边值问题，最普遍的自
10、变量是“位移”。例如，描 述梁弯曲变形的微分方程，边界条件总是规定在 梁的两端。n8.2.1 常微分方程的解析解n在MATLAB中，用大写的字母D表示导数，Dy表示一 阶导数 ；D2y表示二阶导 ；表示微分方程：n在MATLAB中，由函数dsolve进行常微分方程（组 ）的求解问题，其调用格式是：ny=dsolve(e1,e2,.,en,c1,c2,., cn,v1,v2,.,vn)n它用于求解常微分方程组e1,e2,.,en在初值条 件c1,c2,.,cn下的特解。如果不给出初值条件 ，则求出解常微分方程组的通解。v1,v2,.,vn 是求解变量。n例9-1 求常微分方程 的通解n在此问题中
11、，微分方程的MATLAB表达式为：n Dy-y*(1-y2)=0n没有给出初值条件，只能求出解常微分方程组的 通解，调用求解微分方程组的语句为：ny=dsolve(Dy-y*(1-y2)=0,x)n运算结果：ny = 1/(1+exp(-2*x)*C1)(1/2)n -1/(1+exp(-2*x)*C1)(1/2)n即该微分方程的解析解为 n求常微分方程组的通解。n没有给出初值条件，只能求出解常微分方程组的 通解，调用求解微分方程组的语句为：nx,y=dsolve(Dx+y=exp(z),Dy- =exp(2*z),z)n运算结果：nx =ncos(z)*C2-sin(z)*C1-1/5*ex
12、p(2*z)+1/2*exp(z)ny =nsin(z)*C2+cos(z)*C1+2/5*exp(2*z)+1/2*exp(z)n其中，C1与C2是积分常数。即该微分方程的解析 解为 n求常微分方程组当 和 条件下的特解。n在此问题中，两个微分方程的MATLAB表达式为：ne1：Dx+2*x-Dy=10*cos(t)ne2：Dx+Dy+2*y=4*exp(-2*t)n初值条件表达式为：n C1：x(0)=2n C2：y(0)=0 n系统函数dsolve默认自变量为t，所以调用求解微 分方程组的语句为：nx,y=dsolve(Dx+2*x-Dy=10*cos(t),nDx+Dy+2*y=4*e
13、xp(-2*t),x(0)=2,y(0)=0)n运算结果：nx = -2*exp(-t)*sin(t)+4*cos(t)+3*sin(t)- 2*exp(-2*t)ny = 2*exp(-t)*cos(t)+sin(t)-2*cos(t)n即该微分方程的解析解为：n8.2.2 常微分方程的数值解nMATLAB中用于求解一阶微分方程组初值问题数值解 的常用函数是ode45，它是基于四阶五级Runge- Kutta变步长算法，采用误差作为监测指标的闭环 控制，既保证有较高的运算速度和中等计算精度， 也保证预期的数值稳定性，是大部分场合求解非刚 性常微分方程数值解的首选方法。求解非刚性常微 分方程数
14、值解的常用函数还有ode23，它是基于二 阶三级Runge-Kutta的单步算法，计算精度较低。n函数ode45的调用格式是：n t,x=ode45(fun,t0,tf,x0)n其中，fun是描述微分方程的符号函数或是用 inline()函数描述的微分方程；t0表示起始时刻， tf表示终止时刻；x0是初值问题的初始状态变量。n例9-2 试求洛伦茨（lorenz）系统动态模型状态 方程的数值解。n初值 ， ；n终止时刻 。n其中，系数n% 描述洛伦茨（lorenz）系统动态模型状态方程3 维非线性常微分方程组的函数文件nfunction xdot=lorenzeq(t,x,flag,beta,r
15、ho,sigma) n% 输入参数中：beta,rho,sigma是方程组有关项 的系数；变量flag用于占位nxdot=-beta*x(1)+x(2)*x(3);- rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)- x(3);n% 求解洛伦茨（lorenz）系统3维非线性常微分方 程组ntf=100; x0=0;0;1e-10; % tf是终止时刻nbeta=8/3;rho=10;sigma=28; % 给各系数赋值nt,x=ode45(lorenzeq,0,tf,x0,beta,rho,sigma);nsubplot(1,2,1);nplot(t,x) % 绘制2维线图ngrid onntitle(bf 状态变量的时间
