代数方程和常微分方程求解.ppt

第8章 代数方程和常微分方程求解,代数方程是未知数和常数进行有限次代数运算所组成的方程，它包括有理方程和无理方程。代数方程 的解称为 的根或零点，其求解一般是通过代数几何来进行。微分方程是含有一个或是多个导数的方程。只有一个自变量及其导数的微分方程称为常微分方程；包含两个以上自变量及其偏导数的微分方程称为偏微分方程。工程上许多物理规律，设计过程的模拟和评价，凡是涉及质量和能量运动设计分析的问题，都最终归结到微分方程。,8.1 代数方程求解,8.1.1 代数方程图解法符号绘图函数fplot()和ezplot()也可以用于图解法求代数方程的根，它适用于求解维数较少的一维方程或二维方程组。对于一维方程图解，其解就是函数曲线与x轴交点所对应的变量数值。如果有多个交点，则表示该方程有多个解；如果没有交点，则表示该方程没有解。例如，在例5-3使用符号绘图函数绘制代数方程的图形（图5-3左图）中可见，函数在区间-5,5内与x轴有3个交点，因此该代数方程该区间内有3个实根。,对于二维方程组图解，其解就是两条函数曲线的交点所对应的坐标数值。如果只有1个交点（或切点），则表示该方程组有1个解；如果有2个交点，则表示该方程组有2个解；如果没有交点，则表示该方程没有解。例8-1 用图解法求解二维联立方程。a=-2;b=2;%定义横轴区间ezplot(x2+y2-1.69,a,b);axis equal;%控制坐标轴比例相等hold on;grid on;ezplot(2.4*x3-y+1.5,a,b);line(a,b,0,0);line(0,0,b,a);xlabel(bf x);ylabel(bf y);title(bf 二维代数方程组的图解法),gtext(bf f_1=x2+y2-1.32);gtext(bf f_2=2.4x3-y+1.5);,8.1.2 代数方程的解析解求非线性方程或方程组解析解的函数调用格式：X=solve(fun,x)其中，fun是符号方程的函数表达式，x是自变量，X是解析解。应当指出，函数solve(fun,x)也可以用于求线性方组的解析解。例8-2 求非线性解方程组解析解,%二维非线性方程组的解析解syms a b x y;f1=x2-x*y-a;f2=y2-x*y+b;disp(二维非线性方程组的解析解：)X,Y=solve(f1,f2,x,y)M文件运行结果：二维非线性方程组的解析解：x=a/(a-b)(1/2)-a/(a-b)(1/2)Y=1/(a-b)(1/2)*b-1/(a-b)(1/2)*b,8.1.3 线性方程组的数值解最简便方法是使用矩阵左除或是矩阵求逆的方法，求解线性方程组AX=b。X=Ab X=inv(A)*b其中，A是方程组的系数矩阵，b是常数向量，X是解析解。例8-3 求线性方程组的数值解,%线性方程组的数值解AA=1,1,1,1;1,2,-1,4;2,-3,-1,-5;3,1,2,11;bb=5;-2;-2;0;%线性方程组常数向量disp(采用矩阵左除求出线性方程组的解：)xx=AAbbdisp(采用矩阵求逆求出线性方程组的解：)zx=inv(AA)*bbdisp(计算残量：)r=AA*zx-bbdisp(计算残量的模：)R=norm(r),M文件运行结果：采用矩阵左除或矩阵求逆求出线性方程组的解：xx(zx)=1.0000 2.0000 3.0000-1.0000计算残量：r=1.0e-014*0.0888 0.2220-0.4441 0.1776计算残量的模：R=5.3475e-015,8.1.4 非线性方程的数值解1、一维非线性方程对于一维非线性方程求解，可以看作是单变量的极小化问题，通过不断缩小搜索区间来逼近一维问题的真解。因此，可以使用一维非线性方程优化解函数来求解。其调用格式是：x,fx,flag=fzero(fun,x0)其中，输入参数中：fun是非线性方程的函数表达式；x0是根的初值；输出参数中：x是非线性方程的数值解；fx是数值解的函数值；返回参数flag0时，表示求解成功，否则求解出现问题。函数fzero所使用的算法为二分法、secant法和逆二次插值法的组合。,例8-4 求解一维非线性方程%求解单变量x非线性方程x0=0.1;%解的初值xz,fz,flag=fzero(atan(x)+exp(x),x0);disp(求解成功性判断参数：),flagdisp(非线性方程的解：),xzdisp(非线性方程解的函数值：),fzM文件运行结果：求解成功性判断参数：flag=1非线性方程的解：xz=-0.6066非线性方程解的函数值：fz=-1.1102e-016,2、多维非线性方程组对于多维非线性方程组使用多维非线性方程组优化解函数求解的格式：x,fval,flag=fsolve(fun,x0)其中，输入参数中fun是非线性方程组的向量函数表达式；x0是根的初值；输出参数中x是非线性方程（组）的数值解；fval是数值解的函数值；返回参数flag0时，表示求解成功.函数fsolve的作用是从根的初值x0开始，以逐渐减少误差的算法，搜索出满足多维非线性方程组fun的实根x和对应的函数值fval。,例8-5 求解二维非线性方程组x0=2;3;%解的初值%定义非线性方程组表达式f和向量xfun=inline(2*x(1)-x(2)2-exp(-x(1);-x(1)3+x(1)*x(2)-exp(-x(2),x)xn,fval,flag=fsolve(fun,x0);disp(求解成功性判断参数：),flagif flag0 disp(方程组的解成功)elseif flag=0 disp(方程组的解不成功)enddisp(非线性方程组的解：),xndisp(非线性方程组解的函数值：),fval,M文件运行结果：求解成功性判断参数：flag=1方程组的解成功非线性方程组的解：xn=0.9978 1.2755非线性方程组解的函数值：fval=1.0e-006*-0.1945-0.3372两个非线性方程解的函数值非常接近于0，说明其解的误差很小，且判断参数flag0，解成功。,8.2 常微分方程求解,求解微分方程必须事先对自变量的某些值规定出函数或是导数的值。若在自变量为零的点上，给出初始条件，称为初值问题，最普遍的自变量是“时间”。例如，弹性系统的自由振动，若以时间为零来限定位移和速度，这是一个初值问题。若在自变量为非零的点上，给出边界条件，称为边值问题，最普遍的自变量是“位移”。例如，描述梁弯曲变形的微分方程，边界条件总是规定在梁的两端。,8.2.1 常微分方程的解析解在MATLAB中，用大写的字母D表示导数，Dy表示一阶导数；D2y表示二阶导；表示微分方程：在MATLAB中，由函数dsolve进行常微分方程（组）的求解问题，其调用格式是：y=dsolve(e1,e2,.,en,c1,c2,.,cn,v1,v2,.,vn)它用于求解常微分方程组e1,e2,.,en在初值条件c1,c2,.,cn下的特解。如果不给出初值条件，则求出解常微分方程组的通解。v1,v2,.,vn是求解变量。,例9-1 求常微分方程 的通解在此问题中，微分方程的MATLAB表达式为：Dy-y*(1-y2)=0没有给出初值条件，只能求出解常微分方程组的通解，调用求解微分方程组的语句为：y=dsolve(Dy-y*(1-y2)=0,x)运算结果：y=1/(1+exp(-2*x)*C1)(1/2)-1/(1+exp(-2*x)*C1)(1/2)即该微分方程的解析解为,求常微分方程组的通解。没有给出初值条件，只能求出解常微分方程组的通解，调用求解微分方程组的语句为：x,y=dsolve(Dx+y=exp(z),Dy-=exp(2*z),z)运算结果：x=cos(z)*C2-sin(z)*C1-1/5*exp(2*z)+1/2*exp(z)y=sin(z)*C2+cos(z)*C1+2/5*exp(2*z)+1/2*exp(z)其中，C1与C2是积分常数。即该微分方程的解析解为,求常微分方程组 当 和 条件下的特解。在此问题中，两个微分方程的MATLAB表达式为：e1：Dx+2*x-Dy=10*cos(t)e2：Dx+Dy+2*y=4*exp(-2*t)初值条件表达式为：C1：x(0)=2 C2：y(0)=0,系统函数dsolve默认自变量为t，所以调用求解微分方程组的语句为：x,y=dsolve(Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t),x(0)=2,y(0)=0)运算结果：x=-2*exp(-t)*sin(t)+4*cos(t)+3*sin(t)-2*exp(-2*t)y=2*exp(-t)*cos(t)+sin(t)-2*cos(t)即该微分方程的解析解为：,8.2.2 常微分方程的数值解MATLAB中用于求解一阶微分方程组初值问题数值解的常用函数是ode45，它是基于四阶五级Runge-Kutta变步长算法，采用误差作为监测指标的闭环控制，既保证有较高的运算速度和中等计算精度，也保证预期的数值稳定性，是大部分场合求解非刚性常微分方程数值解的首选方法。求解非刚性常微分方程数值解的常用函数还有ode23，它是基于二阶三级Runge-Kutta的单步算法，计算精度较低。函数ode45的调用格式是：t,x=ode45(fun,t0,tf,x0)其中，fun是描述微分方程的符号函数或是用inline()函数描述的微分方程；t0表示起始时刻，tf表示终止时刻；x0是初值问题的初始状态变量。,例9-2 试求洛伦茨（lorenz）系统动态模型状态方程的数值解。初值，；终止时刻。其中，系数,%描述洛伦茨（lorenz）系统动态模型状态方程3维非线性常微分方程组的函数文件function xdot=lorenzeq(t,x,flag,beta,rho,sigma)%输入参数中：beta,rho,sigma是方程组有关项的系数；变量flag用于占位xdot=-beta*x(1)+x(2)*x(3);-rho*x(2)+rho*x(3);-x(1)*x(2)+sigma*x(2)-x(3);%求解洛伦茨（lorenz）系统3维非线性常微分方程组tf=100;x0=0;0;1e-10;%tf是终止时刻,beta=8/3;rho=10;sigma=28;%给各系数赋值t,x=ode45(lorenzeq,0,tf,x0,beta,rho,sigma);subplot(1,2,1);plot(t,x)%绘制2维线图grid ontitle(bf 状态变量的时间响应图)subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3);%绘制3维线图axis(10 42-20 20-20 25);title(bf 状态变量的相空间三维图),8.2.3 高阶微分方程的降阶求解n阶微分方程一般形式：已知输出变量各阶导数初值选择一组状态变量可以把n阶微分方程化成等价的一阶微分方程组这是一个含n个未知函数的一阶微分方程组,初值,例9-3 试求二阶微分方程 的数值解。状态变量初值 和，起始时刻，终止时刻。该二阶微分方程可以通过下面的变换，得到等价的一阶方程组：编制如下函数m文件：function dy=weifen(t,x)dy=zeros(2,1);dy(1)=x(2);%x(2)=zdy(2)=sin(t)-x(2)-x(1);%x(1)=y,然后用函数ode45求解微分方程：t,y=ode45(weifen,0,20,0,6);plot(t,y(:,1);hold on%绘制y的图像plot(t,y(:,2),r-);%绘制y的图像xlabel(bf t);ylabel(bf y);title(bf 二阶微分方程的数值解);legend(实线：y的图像,虚线：Dy的图像)运算结果如图8-3所示。,例9-4 试求二阶微分方程的数值解：该二阶微分方程可以通过下面的变换，得到等价的一阶方程组：已知状态变量的初值起始时刻，终止时刻。%1-计算二阶微分方程的解析解syms ty=dsolve(D2y=-3*cos(2*t)+2*sin(t)+t-3.5,y(0)=0,Dy(0)=0,t),ezplot(y,0 10);%绘制解析解xlabel(bf t);ylabel(bf y);title(bf 二阶微分方程的解析解和数值解)gtext(bf y=(3cos2t)/4-2sint+t3/6-7t2/4+2t-3/4)hold on%2-计算二阶微分方程的数值解fun=inline(x(2);-3*cos(2*t)+2*sin(t)+t-3.5,t,x);t x1=ode45(fun,0,10,0,0);plot(t,x1(:,1),r*);%绘制数值解grid onlegend(曲线：解析解,星号：数值解)运算结果如图8-4所示。,8.2.4 刚性微分方程的数值解刚性方程是指其Jacobian矩阵的特征值相差十分悬殊。在解的性态上表现为，其中一些解变化缓慢，另一些变化快，且相差较悬殊，这类方程常常称为刚性方程，又称为Stiff方程。刚性方程和非刚性方程对解法中步长选择的要求不同。刚性方程一般不适合由ode45这类函数求解，而应该采用ode15s等函数。如果不能分辨是否是刚性方程，先试用函数ode45，再用函数ode15s。求解刚性微分方程数值解常用函数有：ode15s，它是基于Gears反向数值微分的多步算法，具有中等计算精度；ode23s，它是基于Rosebrock的单步算法，计算精度较低。函数ode15s和ode23s的调用格式与函数ode45相同,例9-5 求刚性微分方程组的数值解已知：初始条件 和，起始时刻 和终止时刻。使用函数dsolve可求出解析解为：可见，当时，两个分量均趋于0。下降极快，而下降很慢。,编制描述刚性微分方程组的M文件：%刚性微分方程组的函数文件function fun=stiff(t,y)fun=-0.01*y(1)-99.99*y(2);-100*y(2);调用函数ode15s求解：t,y=ode15s(stiff,0,400,2,1);tstep=length(t);%总步数 minh=min(diff(t);%最小步长 maxh=max(diff(t);%最大步长 plot(t,y);xlabel(bf t);ylabel(bf y);title(bf 刚性微分方程组的数值解)gtext(bf y_1);gtext(bf y_2);,8.2.5 微分代数方程的数值解微分代数方程（differential algebraion equation，DAE）是指微分方程中某些变量受到某些代数方程条件的约束，它的一般形式是：其中，描述 的方法与普通的微分方程相同；为奇异矩阵，由M文件的求解选项变量Mass说明。在无法严格证明是奇异矩阵时，可以先假设它是奇异矩阵，如果在MATLAB求解过程中没有出现 是非奇异矩阵的警告信息，则求解是有效的，否则无效。,例9-6 求微分代数方程组的数值解。已知：初始条件 和，起始时刻 和终止时刻。该方程组的第3个是代数方程，给定3个状态变量之间的约束关系。用矩阵形式表示该微分代数方程：,编写描述微分代数方程的函数文件：function dx=DAE(t,x)dx=-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);x(1)+x(2)+x(3)-1;%Options是用命令odeset设置的可选积分参数，将常系数矩阵M赋给质量矩阵Mass。M=1,0,0;0,1,0;0,0,0;x0=0.8;0.1;0.1;options=odeset;options.Mass=M;t,x=ode15s(DAE,0,20,x0,options);plot(t,x(:,1),r,t,x(:,2),b,t,x(:,3),g);title(微分代数方程的数值解)gtext(bf x_1(t);gtext(bf x_2(t);gtext(bf x_3(t);xlabel(bf t);ylabel(bf x);grid on,