非线性方程和常微分方程的解法

1、实验8 非线性方程和常微分方程的解法一、实验目的学会用MATLAB软件求解非线性方程和常微分方程。二、实验内容与要求 1. 非线性方程的整值解（1）最小二乘法格式：fsolve（fun,）%求方程fun=0在估计值x附近的近似解。例1.72 求方程的解。>>fc=inline(x-exp(-x);>>xl=fsolve(fc,0) xl=0.5671问题1.28：求解方程,观察知有多解，如何求之？先用命令fplot（5*x2*sin(x)-exp(-x),0,0,10）作图1.13，注意5*x2*sin(x)-exp(-x),0中的“，0”是作y=0直线，即x轴。方程在

2、0，10区间从图中可看出有4个解，分别在0，3，6，9附近，所以用命令：>>fun=inline(5*x2*sin(x)-exp(-x);>>fsolve(fun,0,3,6,9,le-6) 得出结果：ans= 0.5918 3.1407 6.2832 9.4248【例 1.73】求解方程组x-0.7sinx-0.2cosy=0先编制函数文件 fu.m:function y=fu (x)y (1)=x (1) - 0.7*sin ( x (1) ) - 0.2*cos( x (2) ) ;y (2)=x (1) - 0.7*cos (x (1) ) + 0.2*sin

3、(x (2) );y = y(1), y(2) ;在命令窗口调用fu运算：>>x1 =fsolve( fu, 0.5,0.5 )x1 =0.5265 0.5079(2) 零点法格式： %求函数在附近的零点。说明：估计值若为标量时，则在附近查找零点，=向量时，则首先要求函数值 ，然后将严格在区间内零点，若找不到，系统将给出提示。【例 1.74】 求函数 的零点。>> fn =inline('sin (x2) / x+x* exp (x) - 4' );>>x=fzero ( fn, 1,2 ) %这里的fn不要加单引号x = 1.0748注意：

4、方程解的估计值 可用fplot作图看出；用function建立函数文件fn，求解调用时fn两边要加单引号，而用inline时fn两边不要加单引号；这两种方法也可解线性方程组。 代数方程的符号解格式： %求解方程，输入参量可是符号表达式或字符表达式。 %对中指定的变量求解方程。 %求解方程组。 %对方程组中指定的n个变量加求解【例1.75】 >>solve ( 'a*x2 + b*x + c')>>solve( 'a*x2 + b*x + c', 'b')>>x,y =solve ( 'x+y = 1&#

5、39;, 'x -11*y= 5')>>a,u,v =solve ('a*u2 + v2', 'u- v =1','a2 -5*a +6')计算结果为：ans = 1/2/a*(-b +(b2-4*a*c)(1/2) ) 1/2/a*(-b -(b2 -4*a*c)(1/2) )ans =-(a*x2+c) /xx = 4/3y =-1/3a =2233u = 1/3+1/3*i*2(1/2) 1/3-1/3*i*2(1/2) 1/4+1/4*i*3(1/2) 1/4-1/4*i*3(1/2) v = -2/3+1/3

6、*i*2(1/2) -2/3-1/3*i*2(1/2) -3/4+1/4*i*3(1/2) -3/4-1/4*i*3(1/2)注意：对于单个的方程或方程组，若不存在符号解，则返回方程（组）的数值解。 问题1.29：用符号法求解问题1.28中的方程，结果不对，所以要验根，多用几种方法相互验证，用符号法解方程，解的表达式不易懂，怎么办？ x =solve('3*x2-exp(x)')x = -2*lambertw (-1/6*3(1/2) ) -2*lambertw (-1,-1/6*3(1/2) ) -2*lambertw (1/6*3(1/2) ) 再用命令： >>

7、vpa (x,3)ans = .912 3.72 -.4603.常微分方程数值解法格式： %区间上，用初始条件求解显示微分方程说明：solver为命令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb,之一。 odefun为显示常微分方程。 tspan 积分区间 （即求解区间） 的向量。要获得问题在其他指的时间点上的解，则令 （要求是单调的）。包含初始条件的向量。求解具体ODE的基本过程如下所示。 根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。F（y, y, y, y(n),t）=0y(0)=y0, y(0)=y1, ,y(n-

8、1)(0)=yn-1而y=y; y(1);y(2); ;y(m-1),n 与m可以不等。 运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2), ,y2=y,y1=y,把高阶（大于2阶）的方程（组）写成一阶微分方程组： 根据与的结果，编写能计算导数的M函数文件odefile. 将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。 因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver,具体如表1.11所示

9、。 表1.11 不同求解器Solver的特点求解器SolverODE类型特点说明ode45非刚性一步算法，4，5阶Runge-Kutta方程，累计截断误差达（x）3大部分场合的首选算法ode23非刚性一步算法，2，3阶Runge-Kutta方程,累计截断误差达（x）3使用于精度较低的情形ode113非刚性多步法，Adams算法，高低精度均可到10-310-6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gears反向数值微分，精度中等若ode45失效时，可尝试使用ode23s刚性一步法，阶Rosebrock算法，低精度当精度较低时，计算时间比ode1

10、5短【例1.76】 求解微分方程y=2y+2x2+2x, 0x0.5, y(0)=1.>> fun=inline('-2*y+2*x2+2*x','x','y');>> x,y=ode23(fun,0,0.5,1);>> x'ans =Columns 1 through 6 0 0.0400 0.0900 0.1400 0.1900 0.2400Columns 7 through 12 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y'ans

11、 =Columns 1 through 7 1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440Columns 8 through 12 0.6222 0.6105 0.6084 0.6154 0.6179【例】求解描述振荡器的经典的Ver der Pol 微分方程分析：令编写函数文件vdp,m: function fy=vdp(t,x)fy=x(2);7*(1-x(1)2)*x(2)-x(1);再在命令窗口中执行：>> Y0=1;0;>> t,x=ode45('vdp',0,40,Y0);>> y

12、=x(:,1);dy=x(:,2);>> plot(t,y,t,dy)图形结果如图.所示. 图.例.图形结果常微分方程的符号解格式：r=dsolve().说明：对给定的常微分方程（组）中指定的符号自变量v，与给定的边界条件和初始条件求符号解（解析解）r若没有指定变量v,则缺省变量为t;在微分方程（组）的表达式eq中，大写字母D表示对自变量（设为x）的微分算子:微分算子D后面的字母则表示为因变量，即待求解的未知函数；初始和边界条件由字符串表示:,等等，分别表示若边界条件少于方程组的阶数，则返回的结果r中会出现任意的常数若该命令找不到解析解，则返回一警告信息，同时返回一空的sym对象，这时，用户可以用命令ode23或ode45求方程组的数值解例1.78>>D1=dsolve(D2y=Dy+exp(x)>> D2=dsolve(Dy)2+y2=1,s)>>D3=dsolve(Dy=a*y,y(0)=b)%带一个初始条件>> D4=dsolve(D2y=-a2*y,y(0)=1,Dy(pi/a)=0) %带两个初始条件>>x,y=