《数学实验第三次讲稿.ppt》由会员分享，可在线阅读，更多相关《数学实验第三次讲稿.ppt（43页珍藏版）》请在三一办公上搜索。
1、数学实验第三讲,方程求解,2023/10/14,1,1)掌握方程求解的三种解法：解析 法、数值解法以及图形表示解的方法;2)学会使用MATLAB软件求解析解、数值解 和图形解；3）通过范例学习怎样建立方程模型和 分析问题的思想。,实验目的,2023/10/14,2,一，建立方程,例子【问题背景】一段时间,美国原子能委员会是按以下方式处理浓缩放射性废物的.他们将废物装入密封性能很好的圆桶中,然后扔到水深300英尺的海里.这种做法是否会造成放射性污染,很自然地引起了生态学家及社会各界的关注.原子能委员会一再保证,圆桶非常坚固,决不会破漏,这种做法是绝对安全的.然而一些工程师们却对此表示怀疑,他们认
2、为圆桶在海底相撞时有可能发生破裂.由此双方展开了一场笔墨官司.究竟谁的意见正确呢?只能让事实说话了!,2023/10/14,3,数学建模的一般步骤,2023/10/14,4,二，方程求解,1，解析方法2，图形放大法3，迭代方法4，区间方法,2023/10/14,5,1，方程求解之解析方法,主要针对一些比较简单的方程以及方程组，比如多项式方程等。同学们以前对方程的求解也是针对这样一些方程进行的。该方法的优点是可以利用纸笔得到简单有效并且精确的解；缺点是可以求解的方程数量太少。Matlab和Maple提供了求方程解析解的函数，可以说对数学演算提供了不少方便。,2023/10/14,6,2，方程求解
3、之图形放大法,图形最大的有点就是直观，试想如果我们有了函数准确的图形，那么曲线和x轴的交点就是我们要求的方程的解。因此我们可以利用图形工具得到方程的解。当然，计算机上的图形不可能等同于函数的真实图形，因为计算机上的图形是曲线上部分点的轨迹而不是全部，因此通过图形不可能得到方程的精确解，甚至它只是一个比较粗略的解，当然，通过对图形的放大可以得到更精确一些的解。同时，这种方法也不适应大量的数据处理。,2023/10/14,7,方程 f(x)=0 1）建立坐标系，画曲线f(x)；2）观察曲线f(x)与x轴相交的交点；3）将其中一个交点进行局部放大；4）该交点的横坐标值就是方程的根。,2图形放大法,2
4、023/10/14,8,例:求方程 x5+2x2+4=0 的一个根.,画方程曲线图（tuxfd.m）x=-3:0.01:3;y=x.5+2*x.2+4;y1=0*x;plot(x,y,x,y1),由此判断：方程的一个根在区间-2，2内，因此将区间-3，3缩小至-2，2，再观察！,该方程有几个根？欲寻找其中一个实根，并且达到一定的精度。,2图形放大法,2023/10/14,9,逐次缩小区间，观察一个根在-1.55-1.5之间。,2图形放大法,2023/10/14,10,3,方程求解之迭代法,迭代法的理论以及方法的出现，对方程求解有着里程碑式的意义。其基本思想如下:,需要求解的方程：f(x)=0（
5、1）经过某种变形得：x=j(x)（2）从而求解方程（1）转化成为求解（2）得不动点。（满足条件x*=j(x*)的点x*称为不动点）为得到方程的不动点，可以构造迭代过程如下：xn+1=j(xn)，n=0，1，x0 定义为迭代初值。,2023/10/14,11,解：第一步 构造迭代函数：x=j(x),例：用迭代方法求解方程 x3-x2-x-1=0。,2023/10/14,12,第二步 迭代设定初值 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
6、/3);%j2(y)z=1+1/z+1/z2;%j3(z)endX，y，z,2023/10/14,13,j1(x)的迭代是失败的（迭代不收敛）。,精确解：x=1.8393,计算结果,2023/10/14,14,迭代函数j2(x)和j3(x)的选取是成功的。精确解为 x=1.8393。并且选取函数j2(x)、j3(x)其收敛速度不一致，前者的速度快些！,结论,1、当遇到迭代不收敛时有什么解决办法？2、如何提高收敛速度？,对于给定的方程 f(x)=0,有多种方式将它改写成等价的形式 x=j(x)。但重要的是如何改写使得序列收敛？,返 回,2023/10/14,15,当今最流行的迭代法是牛顿法以及由
7、此改进的一些方法，比如拟牛顿法等。其基本的思想就是构造迭代格式是利用函数的导数，这类方法有收敛速度快，稳定性好等特点。对低维和高维情况都适合，也是当今一些软件均采用的方法。当然，因为需要函数的导数信息，所以自然对不可微的问题受到制约。该方法的迭代格式为：xk+1=xk-f(xk)/f(xk),2023/10/14,16,例:利用牛顿法求方程x3-x2-x-1=0的根.方法:第一步,给出函数的导函数3x2-2x-1;第二步,给出函数的迭代格式:xk+1=xk-f(xk)/f(xk);设置一定的精度要求,达到即终止.定义函数m文件:function ff=mynewton(x)ff=(x3 x2-
8、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);end,2023/10/14,17,三，解方程函数格式及例子,Matlab对方程的求解提供了以下的一些函数：（1）多项式求根；（2）线性方程组求解；（3）一般的非线性方程（组）求解：,2023/10/14,18,输出:-1.2131-0.9017+0.5753i-0.9017-0.5753i-0.2694+0.9406i-0.2694-0.9406i 0.4168+0.8419i 0.4168-0.8
9、419i 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),roots()语句的用法,2023/10/14,19,4、线性方程组：AX=b 其中A是mn阶矩阵，b是m维向量。,x=A b or x=inv(A)*b,特点：只能求出一个特解。,2023/10/14,20,solution to the following linear system of equations:You can formulate and solve the problem asA=3 11-2;1
10、1-2;1-1 1;b=7;4;19;x=Abx=13.2188-2.3438 3.4375,2023/10/14,21,例:AX=b,解:输入:A=1 2 3;4 5 6;7 8 9;b=6;14;-3;x1=Ab,x2=inv(A)*b输出:警告:系统的秩不足.解不唯一.,1、题中rank(A)=rank(A|b)=23,该方程组有无穷解。2、输出结果是否一致？3、如何求方程组的全部解？,思考,Ab 和inv()语句的用法,返 回,2023/10/14,22,函数fzero,格式：fzero(函数名,初值或区间)Example 1.Calculate by finding the zero
11、 of the sine function near 3.x=fzero(sin,3)x=3.1416Example 2.To find the zero of cosine between 1 and 2 x=fzero(cos,1 2)x=1.5708,2023/10/14,23,1、方程(组)，f1(x)=0,fn(x)=0,x=(x1,xn),2、方程(组)，f1(x)=0,fn(x)=0,x=(x1,xn),solve(f1(x),f2(x),fn(x),X=fsolve(fun,X0,options),MATLAB软件直接求解法,2023/10/14,24,输出:1/2/a*(-b
12、+(b2-4*a*c)(1/2)1/2/a*(-b-(b2-4*a*c)(1/2),单变量方程,solve()语句的用法,例1：求解方程 ax2+bx+c=0,输入:x=solve(a*x2+b*x+c)或 solve(a*x2+b*x+c=0),1）符号解,2023/10/14,25,例2：解方程：x3-2x2=x-1解：s=solve(x3-2*x2=x-1)double(s),2）数字解,该方程是否有实根？vpa(s,10),solve()语句的用法,2023/10/14,26,例3 求解方程：tan(x)-sin(x)=0,3）无穷解,输入:solve(tan(x)-sin(x)=0)
13、输出：0(不能给出全部解）(tx1.m),solve()语句的用法,2023/10/14,27,输入：x,y=solve(x2*y2,x-(y/2)-b)输出：x=0，y=-2*b 0，-2*b（符号解）b，0 b，0 v=x,y,多变量方程组,例4,solve()语句的用法,2023/10/14,28,例6：求解方程组,解 输入：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),fsolve()语句的用法,2023/10/14,29,解：1）建立方程组的M-函数文件(nxxf.m)function eq=nxxf(
14、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;,运行程序(test4.m)y=fsolve(nxxf,1,1,1),3）运行结果：Optimization Terminated Successfully y=0.5990 2.3959 2.0050,fsolve()语句的用法,2023/10/14,30,fsolve()函数的第三个输入是options,它是一个结构型数据,可以通过函数optimset()进行设定.当不进行设定时采用缺省设置.fsolve函数还可以有后面的参
15、数设定,这一功能在有的时候非常有用.,fsolve()语句的用法,例如,求解函数sin(ax)-x=0的最小正解.分析:已知的知识我们知道,该方程没有解析形式的解,也就是说很难得到这个解和a之间的具体关系.数值的方法可以对a先取确定的值,这样变成一个一元方程,容易进行求解,当我们变化a时,就得到很多这样的解,通过这样的方法,我们可以得到解和a之间的一些大致关系.,2023/10/14,31,具体方法如下:建立函数m文件:function ff=funpara(x,a)ff=sin(a*x)-x;,fsolve()语句的用法,建立相应的命令m文件:B=zeros(100,1);for a=1:1
16、00 x0=pi/(2*a)+0.01;B(a)=fsolve(funpara,x0,a);endplot(B),2023/10/14,32,2023/10/14,33,例：某物体的边缘呈圆形，通过测量边沿上是一个点的坐标，数据为（见程序）。使用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.13
17、13 2.4713-0.3435-2.3887-2.9927-1.9572 0.3778 3.2455 5.7042 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点
18、，然后通过缩小区间得到解的近似。优点是不需要函数的导数信息，而且只要有0点就一定可以得到；缺点是相对牛顿法等速度较慢。,2023/10/14,34,问题的关键在于圆桶到底能承受多大速度的碰撞?圆桶和海底碰撞时的速度有多大?工程师们进行了大量破坏性的实验,发现圆桶在直线速度为40 ft/s 的冲撞下会发生破裂,剩下的问题就是计算圆桶沉入300 ft 深的海底时,其末速度究竟有多大?,问题分析,引例的分析和求解,2023/10/14,35,1.使用55加仑的圆桶;(1加仑=3.7854升)2.装满放射性废物时的圆桶重量为 W=527.436磅(1 磅=0.4526公斤)3.在海水中圆桶受到的浮力
19、B=470.327磅 4.圆桶下沉时受到海水的阻力 D=C v C 为常数,经测算得:C=0.08.5.建立坐标系,取垂直向下为坐标方向 y,海平面为坐标原点.,问题假设,引例,2023/10/14,36,根据牛顿第二定律,圆桶下沉时应满足微分方程:,建立模型,引例,2023/10/14,37,为了求出圆桶与海底的碰撞速度v(t),需要求出圆桶下沉到海底300英尺时的时间 t,再计算v(t)，要做到这一点是十分困难的.若将速度v 看成是海水深度y 的函数,即,由复合函数的求导法知,建立模型,引例,2023/10/14,38,借助数值方法求出 v(300)=45.1ft/s，显然大于40ft/s
20、。,结论：放射性废物不能随意放入公海！,建立模型,引例,2023/10/14,39,1、求下列方程的根 1）xsin(x)+cos(x)-sin(x)-2x=0,在-1，1上的近似解。2）判定方程 x7+x5+1=0 有几个实根?3)找出函数 f(x)=x3-6x2-2x+12 的所有零点。,2、求解线性方程组：,实 验 内 容,2023/10/14,40,3、求下列方程组的解,使用命令 solve()或 fsolve()。,x0=-5,-5,实 验 内 容,2023/10/14,41,4.炮弹发射视为斜抛运动，已知初始速度为 200 m/s，问要击中水平距离360m、垂直距离160m 的目标
21、，当忽略空气阻力时，发射角应多大？进一步思考：如果要考虑水平方向的阻力，且设阻力与（水平方向）速度成正比，系数为 0.1（1/s），结果又如何？炮弹质量假设为10千克？,实 验 内 容,2023/10/14,42,5.讨论:有这样的电视节目,对一件商品要求观众猜价格,主持人对观众所猜的数目的回答是高或者低(相对于实际价格).现在的问题时,对于价格在1000元之内的商品,寻找一个最好的方法,保证在最少的次数猜中商品价格.试给出竞猜5次得到正确答案的概率位多大.(设商品最小价格单位为元),实 验 内 容,结 束,提示:利用二分法的思想,可以比较快的得到竞猜次数的最大值.一个更有趣的方法是二叉树,其深度就是对应的解.对于平均次数,可以采用模拟的方法,比如产生10000次随机数,看每次猜中的次数,平均就得到平均次数.,2023/10/14,43,
