数学实验第三次讲稿

1、2022-4-131 1) 掌握方程求解的三种解法：解析 法、数值解法以及图形表示解的方法; 2) 学会使用MATLAB软件求解析解、数值解 和图形解； 3） 通过范例学习怎样建立方程模型和 分析问题的思想。实验目的2022-4-132一，引例 例子【问题背景问题背景】一段时间一段时间, 美国原子能委员会是美国原子能委员会是按以下方式处理浓缩放射性废物的按以下方式处理浓缩放射性废物的. 他们将废物装入他们将废物装入密封性能很好的圆桶中密封性能很好的圆桶中, 然后扔到水深然后扔到水深300英尺的海英尺的海里里. 这种做法是否会造成放射性污染这种做法是否会造成放射性污染, 很自然地引起很自然地引起

2、了生态学家及社会各界的关注了生态学家及社会各界的关注. 原子能委员会一再保原子能委员会一再保证证, 圆桶非常坚固圆桶非常坚固, 决不会破漏决不会破漏, 这种做法是绝对安这种做法是绝对安全的全的. 然而一些工程师们却对此表示怀疑然而一些工程师们却对此表示怀疑, 他们认为他们认为圆桶在海底相撞时有可能发生破裂圆桶在海底相撞时有可能发生破裂. 由此双方展开了由此双方展开了一场笔墨官司一场笔墨官司. 究竟谁的意见正确呢究竟谁的意见正确呢? 只能让事实只能让事实说话了说话了!2022-4-133数学建模的一般步骤模型准备模型准备模型假设模型假设模型构成模型构成模型检验模型检验模型分析模型分析模型求解模型

3、求解模型应用模型应用2022-4-134二，方程求解 1，解析方法 2，图形放大法 3，迭代方法 4，区间方法2022-4-1351，方程求解之解析方法 主要针对一些比较简单的方程以及方程组，比如多项式方程等。同学们以前对方程的求解也是针对这样一些方程进行的。该方法的优点是可以利用纸笔得到简单有效并且精确的解；缺点是可以求解的方程数量太少。Matlab和Maple提供了求方程解析解的函数，可以说对数学演算提供了不少方便。2022-4-1362，方程求解之图形放大法 图形最大的有点就是直观，试想如果我们有了函数准确的图形，那么曲线和x轴的交点就是我们要求的方程的解。因此我们可以利用图形工具得到方

4、程的解。当然，计算机上的图形不可能等同于函数的真实图形，因为计算机上的图形是曲线上部分点的轨迹而不是全部，因此通过图形不可能得到方程的精确解，甚至它只是一个比较粗略的解，当然，通过对图形的放大可以得到更精确一些的解。同时，这种方法也不适应大量的数据处理。2022-4-137 方程 f(x)=0 1）建立坐标系，画曲线f(x)； 2）观察曲线f(x)与x轴相交的交点； 3）将其中一个交点进行局部放大； 4）该交点的横坐标值就是方程的根。 2022-4-138例: 求方程 x5 +2x2 + 4 = 0 的一个根.画方程曲线图（画方程曲线图（tuxfd.m）x=-3:0.01:3;y=x.5+2*

5、x.2+4;y1=0*x;plot(x,y,x,y1)-6-4-20246-8000-6000-4000-200002000400060008000由此判断：方程的一个根在区间由此判断：方程的一个根在区间-2-2，22内，因此将区间内，因此将区间-3-3，33缩小至缩小至-2-2，22，再观察！，再观察！ 该方程有几个根？欲寻找其中一个实根，并且达到一定的精度。2022-4-139 逐次缩逐次缩小区间，观小区间，观察一个根在察一个根在-1.55-1.5之之间。间。-202-50050-2-1.5-1-20-10010-2-1.5-20-10010-1.6-1.5-1.4-2-1012022-4

6、-13103,方程求解之迭代法 迭代法的理论以及方法的出现，对方程求解有着里程碑式的意义。其基本思想如下: 需要求解的方程： f (x) = 0 （1） 经过某种变形得：x = j (x) （2） 从而求解方程（1）转化成为求解（2）得不动点。（满足条件x*=j(x*)的点x*称为不动点） 为得到方程的不动点，可以构造迭代过程如下： xn+1 =j (xn)，n =0，1， x0 定义为迭代初值。2022-4-1311解：解： 第一步第一步 构造迭代函数： x=j (x)x(x1x11x)x(1xxx)x(1xxx32232123j+j+j例：例：用迭代方法求解方程 x3 x2 x1 0。20

7、22-4-1312第二步第二步 迭代设定初值 x0=1， xn+1 =j (xn)，n =0，1， 用 MATLAB 编程（died2.m文件)x=1;y=1;z=1;(初始点）for k=1:20 x=x3-x2-1; % j1 (x) y=(y2+y+1)(1/3); % j2 (y) z=1+1/z+1/z2; % j3 (z)endX，y，z2022-4-1313序号j3(x)序号j2(x)j3(x)1j2(x)1.44223.000081.81751.813621.65371.444491.83851.855431.75322.1716101.83891.829441.79951.6

8、725111.83911.845451.82091.9554121.83921.835561.83081.7730131.83921.841671.83541.8822j1(x)的迭代是失败的（迭代不收敛 ）。精确解：x=1.8393计算结果计算结果2022-4-1314 迭代函数j2(x)和j3(x)的选取是成功的。精确解为 x=1.8393。 并且选取函数j2(x)、j3(x)其收敛速度不一致，前者的速度快些！结论1、当遇到迭代不收敛时有什么解决办法？2、如何提高收敛速度？ 对于给定的方程 f(x) = 0, 有多种方式将它改写成等价的形式 x = j(x)。但重要的是如何改写使得序列如何

9、改写使得序列收敛？收敛？2022-4-1315 当今最流行的迭代法是牛顿法以及由此改进的一些方法，比如拟牛顿法等。其基本的思想就是构造迭代格式是利用函数的导数，这类方法有收敛速度快，稳定性好等特点。对低维和高维情况都适合，也是当今一些软件均采用的方法。当然，因为需要函数的导数信息，所以自然对不可微的问题受到制约。该方法的迭代格式为： xk+1=xk-f(xk)/f(xk)2022-4-1316例:利用牛顿法求方程x3 -x2 -x-1 = 0的根.方法:第一步,给出函数的导函数3x2-2x-1; 第二步,给出函数的迭代格式: xk+1=xk-f(xk)/f(xk); 设置一定的精度要求,达到即

10、终止.定义函数m文件: function ff=mynewton(x) ff=(x3 x2 -x-1 )/(3*x2-2*x-1);定义命令m文件: x0=1; x1=x0-mynewton (x0); while abs(x1-x0)0.0001 x0=x1; x1=x0-mynewton (x0); end2022-4-1317三，解方程函数格式及例子 Matlab对方程的求解提供了以下的一些函数： （1）多项式求根； （2） 线性方程组求解； （3）一般的非线性方程（组）求解：2022-4-1318-1.3-1.2-1.1-1-1.5-1-0.500.51x0输出: -1.2131 -0

11、.9017 + 0.5753i -0.9017 - 0.5753i -0.2694 + 0.9406i -0.2694 - 0.9406i 0.4168 + 0.8419i 0.4168 - 0.8419i 0.8608 + 0.3344i 0.8608 - 0.3344i例：求解多项式方程 x9+x8+1=0输入: p=1,1,0,0,0,0,0,0,0,1; roots(p)1，roots()语句的用法-10-50510-2-1.5-1-0.500.511.52x 108x02022-4-1319线性方程组：AX = bAX = b 其中A是mn阶矩阵，b是m维向量。 x=A b or x

12、=inv(A)*b特点：只能求出一个特解。2022-4-13202，线性方程组solution to the following linear system of equations: You can formulate and solve the problem as A = 3 11 -2; 1 1 -2; 1 -1 1; b = 7; 4; 19; x = Ab x = 13.2188 -2.3438 3.43752022-4-13213146,987654321Ab例: AX = b, 解: 输入:A=1 2 3; 4 5 6; 7 8 9; b=6; 14; -3; x1=Ab, x

13、2=inv(A)*b输出:警告: 系统的秩不足. 解不唯一.1、题中rank(A)=rank(A|b)=23, 该方程组有无穷解。2、输出结果是否一致？3、如何求方程组的全部解？思考返 回2022-4-13223.1，函数fzero的用法格式： fzero(函数名,初值或区间) Example 1. Calculate by finding the zero of the sine function near 3. x = fzero(sin,3) x = 3.1416 Example 2. To find the zero of cosine between 1 and 2 x = fzer

14、o(cos,1 2) x = 1.57082022-4-13231、方程(组)， f1(x) = 0,fn(x) = 0, x = (x1,xn)2、方程(组)， f1(x) = 0,fn(x) = 0, x = (x1,xn)fun.m function f = fun(x) f(1)= f1(x) ; f(n)= fn(x) 初值1）可以省略。2）options=1,表示输出中间结果。solve(f1(x),f2(x),fn(x) ) X = fsolve (fun, X0, options) MATLAB软件直接求解法2022-4-1324输出: 1/2/a*(-b+(b2-4*a*c)

15、(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2)单变量方程单变量方程0)(xf3.2 solve()语句的用法例1： 求解方程 ax2+bx+c = 0输入: x= solve(a*x2+b*x+c)或 solve(a*x2+b*x+c=0)1）符号解2022-4-1325例2： 解方程： x3-2x2=x-1解： s=solve(x3-2*x2=x-1) double(s)2）数字解该方程是否有实根？vpa(s,10) 2022-4-1326 例3 求解方程： tan(x)-sin(x)=0 3）无穷解输入: solve(tan(x)-sin(x)=0)输出：0 (不能给出全部

16、解） (tx1.m)2022-4-1327输入： x,y=solve(x2*y2, x-(y/2)-b)输出：x = 0 ， y = -2*b 0， -2*b （符号解） b， 0 b， 0 v=x,y 多变量方程组多变量方程组0)(,0)(1xfxfm例4byxyx2022byxyx2122或 b=2solve()语句的用法2022-4-1328+05012307lnsin32zyxzxzyxy例6：求解方程组解 输入：syms x,y,z;x,y,z=solve(sin(x)+y2+log(z)-7=0, 3*x+2y-z3+1=0,x+y+z-5=0,x,y,z)x =0.5990537

17、5664056731520568183824539y =2.3959314023778168490940003756591z =2.0050148409816158357003177860955输出：3.3 fsolve()语句的用法2022-4-1329解：1）建立方程组的M-函数文件(nxxf.m) function eq=nxxf(x) eq(1)=sin(x(1)+x(2)2+log(x(3)-7; eq(2)=3*x(1)+2x(2)-x(3)3+1; eq(3)=x(1)+x(2)+x(3)-5;2) 运行程序(test4.m) y=fsolve(nxxf,1,1,1)3）运行结果

18、：Optimization Terminated Successfully y= 0.5990 2.3959 2.0050fsolve()语句的用法2022-4-1330fsolve()函数的第三个输入是函数的第三个输入是options,它是一个结构型数据它是一个结构型数据,可以可以通过函数通过函数optimset ()进行设定进行设定.当不进行设定时采用缺省设置当不进行设定时采用缺省设置. fsolve函数还可以有后面的参数设定函数还可以有后面的参数设定,这一功能在有的时候非常这一功能在有的时候非常有用有用.fsolve()语句的用法例如,求解函数sin(ax)-x=0的最小正解.分析: 已

19、知的知识我们知道,该方程没有解析形式的解,也就是说很难得到这个解和a之间的具体关系.数值的方法可以对a先取确定的值,这样变成一个一元方程,容易进行求解,当我们变化a时,就得到很多这样的解,通过这样的方法,我们可以得到解和a之间的一些大致关系.2022-4-1331具体方法如下具体方法如下:建立函数建立函数m文件文件:function ff=funpara(x,a)ff=sin(a*x)-x;fsolve()语句的用法建立相应的命令m文件:B=zeros(100,1); for a=1:100 x0=pi/(2*a)+0.01;B(a)=fsolve(funpara,x0,a);endplot(

20、B)2022-4-13322022-4-1333例：例：某物体的边缘呈圆形，通过测量边沿上是一个点的坐标，数据为（见程序）。使用fsolve计算物体边缘的曲线方程。解：假定物体边缘的曲线方程为（x-a）2+(y-b)2=r2为了得到a，b和r的值。将上述十一组数据代人方程，得到（xi-a）2+(yi-b)2=r2（i=1,11）。次数方程个数是11个，未知量个数是3个。使用fsolve求解，其函数m文件为：function y=funcircle(x)A=6.7630 5.1313 2.4713 -0.3435 -2.3887 -2.9927 -1.9572 0.3778 3.2455 5.7

21、042 6.9465;B=23.2879 25.6492 26.7268 26.1668 24.1531 21.3470 18.6699 17.0010 16.8883 18.3688 20.9564;for i=1:11 y(i)=(A(i)-x(1)2+(B(i)-x(2)2-x(3)2;end二分法 对于一个闭区间上的连续函数，我们有一个0点存在定理，利用这个定理，可以对不可微的函数求得函数的0点。 基本思想是通过判断函数在端点处的函数值异号可以确定函数在开区间上至少有一个0点，然后通过缩小区间得到解的近似。 优点是不需要函数的导数信息，而且只要有0点就一定可以得到；缺点是相对牛顿法等速

22、度较慢。2022-4-1334 问题的关键在于圆桶到底能承受多大速度问题的关键在于圆桶到底能承受多大速度的碰撞的碰撞? 圆桶和海底碰撞时的速度有多大圆桶和海底碰撞时的速度有多大? 工程师们进行了大量破坏性的实验工程师们进行了大量破坏性的实验, 发现发现圆桶在直线速度为圆桶在直线速度为40 ft/s 的冲撞下会发生破裂的冲撞下会发生破裂, 剩下的问题就是计算圆桶沉入剩下的问题就是计算圆桶沉入300 ft 深的海底深的海底时时, 其末速度究竟有多大其末速度究竟有多大? 问题分析问题分析引例的分析和求解2022-4-1335 1. 使用使用55加仑的圆桶加仑的圆桶; ( 1加仑加仑 = 3.7854

23、升升 ) 2. 装满放射性废物时的圆桶重量为装满放射性废物时的圆桶重量为 W = 527.436磅磅 (1 磅磅 = 0.4526公斤公斤 ) 3. 在海水中圆桶受到的浮力在海水中圆桶受到的浮力 B = 470.327磅磅 4. 圆桶下沉时受到海水的阻力圆桶下沉时受到海水的阻力 D = C v C 为常数为常数, 经测算得经测算得: C = 0.08. 5. 建立坐标系建立坐标系, 取垂直向下为坐标方向取垂直向下为坐标方向 y , 海平面为坐标原点海平面为坐标原点.y0 问题假设问题假设引例2022-4-1336根据牛顿第二定律根据牛顿第二定律, , 圆桶下沉时应满足微分方程圆桶下沉时应满足微

24、分方程: :)(/(86.713)(lim)1 ()(0)0(:dd,dd22时当速度：容易计算出圆桶的极限其解：初值条件其中阻力）浮力（重力tsftCBWtveCBWtvvBCvWdtdvmvtyCvDgWmDBWtymtWCg327.470B 建立模型建立模型引例2022-4-1337 为了求出圆桶与海底的碰撞速度为了求出圆桶与海底的碰撞速度v(t), 需要求出需要求出圆桶下沉到海底圆桶下沉到海底300英尺时的时间英尺时的时间 t, 再计算再计算v(t)，要要做到这一点是十分困难的做到这一点是十分困难的. 若将速度若将速度v 看成是海水深看成是海水深度度y 的函数的函数, 即即vyvtyyvtvtytyvtvdddddddddd)()(22由复合函数的求导法知由复合函数的求导法知 建立模型建立模型引例2022-4-1338WgyBWCvBWCBWCvyvWgyvCvBWvCvBWyvmvln0)0(,0)0(dddd2积分，得：初值条件：或微分方程变为：借助数值方