《代数方程和常微分方程求解.ppt》由会员分享,可在线阅读,更多相关《代数方程和常微分方程求解.ppt(40页珍藏版)》请在金锄头文库上搜索。
1、第第8 8章章 代数方程和常微分方程求解代数方程和常微分方程求解n代数方程是未知数和常数进行有限次代数运算所组成的方程,它包括有理方程和无理方程。代数方程 的解称为 的根或零点,其求解一般是通过代数几何来进行。n微分方程是含有一个或是多个导数的方程。只有一个自变量及其导数的微分方程称为常微分方程;包含两个以上自变量及其偏导数的微分方程称为偏微分方程。n工程上许多物理规律,设计过程的模拟和评价,凡是涉及质量和能量运动设计分析的问题,都最终归结到微分方程。8.1 8.1 代数方程求解代数方程求解n8.1.1 8.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-18-1 用图解法求解二维联立方程。na=
3、-2;b=2; % 定义横轴区间nezplot(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 8.1.2 代数方程的解析解代数方程的解析解n求非线性方程或方程组解析解的函数调用格式:n X=solv
4、e(fun,x)n其中,fun是符号方程的函数表达式,x是自变量,X是解析解。n应当指出,函数solve(fun,x)也可以用于求线性方组的解析解。n例例8-2 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 8.1.3
5、 线性方程组的数值解线性方程组的数值解n最简便方法是使用矩阵左除或是矩阵求逆的方法,求解线性方程组AX=b。n X= Abn X=inv(A)*bn其中,A是方程组的系数矩阵,b是常数向量,X是解析解。n例例8-3 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
6、( 计算残量的模:)nR=norm(r)nM文件运行结果:n采用矩阵左除或矩阵求逆求出线性方程组的解:nxx (zx)= 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 8.1.4 非线性方程的数值解非线性方程的数值解n1、一维非线性方程n对于一维非线性方程求解,可以看作是单变量的极小化问题,通过不断缩小搜索区间来逼近一维问题的真解。因此,可以使用一维非线性方程优化解函数来求解。其调用格式是:n x,fx,fla
7、g=fzero(fun,x0)n其中,输入参数中:fun是非线性方程的函数表达式;x0是根的初值;n输出参数中:x是非线性方程的数值解;fx是数值解的函数值;返回参数flag0时,表示求解成功,否则求解出现问题。n函数fzero所使用的算法为二分法、secant法和逆二次插值法的组合。n例例8-4 8-4 求解一维非线性方程n% 求解单变量x非线性方程nx0=0.1; % 解的初值nxz,fz,flag=fzero(atan(x)+exp(x),x0);ndisp( 求解成功性判断参数:), flagndisp( 非线性方程的解:),xzndisp( 非线性方程解的函数值:),fznM文件运行
8、结果:n求解成功性判断参数:nflag = 1n非线性方程的解:nxz = -0.6066n非线性方程解的函数值:nfz = -1.1102e-016n2、多维非线性方程组n对于多维非线性方程组n使用多维非线性方程组优化解函数求解的格式:n x,fval,flag=fsolve(fun,x0)n其中,输入参数中fun是非线性方程组的向量函数表达式;x0是根的初值;输出参数中x是非线性方程(组)的数值解;fval是数值解的函数值;返回参数flag0时,表示求解成功.n函数fsolve的作用是从根的初值x0开始,以逐渐减少误差的算法,搜索出满足多维非线性方程组fun的实根x和对应的函数值fval。
9、n例例8-5 8-5 求解二维非线性方程组nx0=2;3; % 解的初值n% 定义非线性方程组表达式f和向量xnfun=inline(2*x(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 8.2 常微分方程求解常微分方程求解n求解微分方程必须事先对自变量的某些值规定出函数或是导数的值。n若在自变量为零的点上,给出初始条件,称为初值问题,最普遍的自变
10、量是“时间”。例如,弹性系统的自由振动,若以时间为零来限定位移和速度,这是一个初值问题。n若在自变量为非零的点上,给出边界条件,称为边值问题,最普遍的自变量是“位移”。例如,描述梁弯曲变形的微分方程,边界条件总是规定在梁的两端。n8.2.1 8.2.1 常微分方程的解析解常微分方程的解析解n在MATLAB中,用大写的字母D表示导数,Dy表示一阶导数 ;D2y表示二阶导 ; 表示微分方程:n在MATLAB中,由函数dsolve进行常微分方程(组)的求解问题,其调用格式是:ny=dsolve(e1,e2,.,en,c1,c2,.,cn,v1,v2,.,vn)n它用于求解常微分方程组e1,e2,.,
11、en在初值条件c1,c2,.,cn下的特解。如果不给出初值条件,则求出解常微分方程组的通解。v1,v2,.,vn是求解变量。n例例9-1 9-1 求常微分方程 的通解n在此问题中,微分方程的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没有给出初值条件,只能求出解常微分方程组的通解,调用求解
12、微分方程组的语句为:nx,y=dsolve(Dx+y=exp(z),Dy-=exp(2*z),z)n运算结果:nx =ncos(z)*C2-sin(z)*C1-1/5*exp(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
13、)=0 n系统函数dsolve默认自变量为t,所以调用求解微分方程组的语句为:nx,y=dsolve(Dx+2*x-Dy=10*cos(t),nDx+Dy+2*y=4*exp(-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 8.2.2 常微分方程的数值解常微分方程的数值解nMATLAB中用于求解一阶微分方程组初值问题数值解的常用函数是ode45,它是基于四阶五级Runge-K
14、utta变步长算法,采用误差作为监测指标的闭环控制,既保证有较高的运算速度和中等计算精度,也保证预期的数值稳定性,是大部分场合求解非刚性常微分方程数值解的首选方法。求解非刚性常微分方程数值解的常用函数还有ode23,它是基于二阶三级Runge-Kutta的单步算法,计算精度较低。n函数ode45的调用格式是:n t,x=ode45(fun,t0,tf,x0)n其中,fun是描述微分方程的符号函数或是用inline()函数描述的微分方程;t0表示起始时刻,tf表示终止时刻;x0是初值问题的初始状态变量。n例例9-29-2 试求洛伦茨(lorenz)系统动态模型状态方程的数值解。n初值 , ;n终
15、止时刻 。n其中,系数n% 描述洛伦茨(lorenz)系统动态模型状态方程3维非线性常微分方程组的函数文件nfunction xdot=lorenzeq(t,x,flag,beta,rho,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;sig
16、ma=28; % 给各系数赋值nt,x=ode45(lorenzeq,0,tf,x0,beta, rho,sigma);nsubplot(1,2,1);nplot(t,x) % 绘制2维线图ngrid onntitle(bf 状态变量的时间响应图)nsubplot(1,2,2);nplot3(x(:,1),x(:,2),x(:,3); % 绘制3维线图naxis(10 42 -20 20 -20 25);ntitle(bf 状态变量的相空间三维图)n8.2.3 8.2.3 高阶微分方程的降阶求解高阶微分方程的降阶求解nn阶微分方程一般形式:n已知输出变量各阶导数初值n选择一组状态变量n可以把n
17、阶微分方程化成等价的一阶微分方程组n这是一个含n个未知函数的一阶微分方程组,初值n例例9-39-3 试求二阶微分方程 的数值解。状态变量初值 和 ,起始时刻 ,终止时刻 。n该二阶微分方程可以通过下面的变换,得到等价的一阶方程组:n编制如下函数m文件:nfunction dy=weifen(t,x)ndy=zeros(2,1);ndy(1)=x(2); % x(2)=zndy(2)=sin(t)-x(2)-x(1); % x(1)=y n然后用函数ode45求解微分方程:t,y=ode45(weifen,0,20,0,6);nplot(t,y(:,1);hold on % 绘制y的图像nplo
18、t(t,y(:,2), r-); % 绘制y的图像nxlabel(bf t);nylabel(bf y);ntitle(bf 二阶微分方程的数值解);nlegend(实线:y的图像,虚线:Dy的图像)n运算结果如图8-3所示。n例例9-49-4 试求二阶微分方程的数值解:n该二阶微分方程可以通过下面的变换,得到等价的一阶方程组:n已知状态变量的初值n起始时刻 ,终止时刻 。n% 1-计算二阶微分方程的解析解nsyms tny=dsolve(D2y=-3*cos(2*t)+2*sin(t)+t-3.5,y(0)=0,Dy(0)=0,t)nezplot(y,0 10); % 绘制解析解nxlabe
19、l(bf t);ylabel(bf y);ntitle(bf 二阶微分方程的解析解和数值解)ngtext(bf y=(3cos2t)/4-2sint+t3/6-7t2/4+2t-3/4)nhold on% 2-计算二阶微分方程的数值解fun=inline(x(2);-3*cos(2*t)+2*sin(t)+t-3.5,t,x); nt x1=ode45(fun,0,10,0,0);nplot(t,x1(:,1),r*); % 绘制数值解ngrid onlegend(曲线:解析解,星号:数值解)n运算结果如图8-4所示。 n8.2.4 8.2.4 刚性微分方程的数值解刚性微分方程的数值解n刚性方
20、程是指其Jacobian矩阵的特征值相差十分悬殊。在解的性态上表现为,其中一些解变化缓慢,另一些变化快,且相差较悬殊,这类方程常常称为刚性方程,又称为Stiff方程。n刚性方程和非刚性方程对解法中步长选择的要求不同。刚性方程一般不适合由ode45这类函数求解,而应该采用ode15s等函数。如果不能分辨是否是刚性方程,先试用函数ode45,再用函数ode15s。求解刚性微分方程数值解常用函数有:node15s,它是基于Gears反向数值微分的多步算法,具有中等计算精度;node23s,它是基于Rosebrock的单步算法,计算精度较低。n函数ode15s和ode23s的调用格式与函数ode45相
21、同 n例例9-59-5 求刚性微分方程组的数值解n已知:初始条件 和 ,起始时刻 和终止时刻 。n使用函数dsolve可求出解析解为:n可见,当时,两个分量均趋于0。下降极快,而下降很慢。n编制描述刚性微分方程组的M文件:n% 刚性微分方程组的函数文件nfunction fun=stiff(t,y)nfun=-0.01*y(1)-99.99*y(2);-100*y(2);n调用函数ode15s求解:nt,y=ode15s(stiff,0,400,2,1); ntstep=length(t); % 总步数 nminh=min(diff(t); % 最小步长 nmaxh=max(diff(t);
22、% 最大步长 nplot(t,y); nxlabel(bf t);ylabel(bf y);ntitle(bf 刚性微分方程组的数值解)ngtext(bf y_1); gtext(bf y_2);n8.2.5 8.2.5 微分代数方程的数值解微分代数方程的数值解n微分代数方程(differential algebraion equation,DAE)是指微分方程中某些变量受到某些代数方程条件的约束,它的一般形式是:n其中,描述 的方法与普通的微分方程相同; 为奇异矩阵,由M文件的求解选项变量Mass说明。n在无法严格证明是奇异矩阵时,可以先假设它是奇异矩阵,如果在MATLAB求解过程中没有出现
23、 是非奇异矩阵的警告信息,则求解是有效的,否则无效。n例例9-69-6 求微分代数方程组的数值解。已知:初始条件 和 ,起始时刻 和终止时刻 。n该方程组的第3个是代数方程,给定3个状态变量之间的约束关系。n用矩阵形式表示该微分代数方程:n编写描述微分代数方程的函数文件:nfunction dx=DAE(t,x)ndx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);n 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);n x(1)+x(2)+x(3)-1;n% Options是用命令odeset设置的可选积分参数,将常系数矩阵M赋给质量矩阵Mass。nM=1,0,0;0,1,0;0,0,0;x0=0.8;0.1;0.1;noptions=odeset;options.Mass=M;nt,x=ode15s(DAE,0,20,x0,options);nplot(t,x(:,1),r,t,x(:,2),b, t,x(:,3),g); title(微分代数方程的数值解)ngtext(bf x_1(t);gtext(bf x_2(t);ngtext(bf x_3(t); xlabel(bf t);ylabel(bf x);grid on