数学实验第三次讲稿市公开课一等奖百校联赛特等奖课件.pptx

/10/101第1页 1)掌握方程求解三种解法：解析 法、数值解法以及图形表示解方法;2)学会使用MATLAB软件求解析解、数值解 和图形解；3）经过范例学习怎样建立方程模型和 分析问题思想。实验目/10/102第2页一，建立方程例子【问题背景问题背景】一段时间一段时间,美国原子能委员会是美国原子能委员会是按以下方式处理浓缩放射性废物按以下方式处理浓缩放射性废物.他们将废物装入他们将废物装入密封性能很好圆桶中密封性能很好圆桶中,然后扔到水深然后扔到水深300英尺海英尺海里里.这种做法是否会造成放射性污染这种做法是否会造成放射性污染,很自然地引发很自然地引发了生态学家及社会各界关注了生态学家及社会各界关注.原子能委员会一再保原子能委员会一再保证证,圆桶非常坚固圆桶非常坚固,决不会破漏决不会破漏,这种做法是绝对安这种做法是绝对安全全.然而一些工程师们却对此表示怀疑然而一些工程师们却对此表示怀疑,他们认为他们认为圆桶在海底相撞时有可能发生破裂圆桶在海底相撞时有可能发生破裂.由此双方展开了由此双方展开了一场笔墨官司一场笔墨官司.终究谁意见正确呢终究谁意见正确呢?只能让事实只能让事实说话了说话了!/10/103第3页数学建模普通步骤模型准备模型准备模型假设模型假设模型组成模型组成模型检验模型检验模型分析模型分析模型求解模型求解模型应用模型应用/10/104第4页二，方程求解1，解析方法2，图形放大法3，迭代方法4，区间方法/10/105第5页1，方程求解之解析方法主要针对一些比较简单方程以及方程组，比如多项式方程等。同学们以前对方程求解也是针对这么一些方程进行。该方法优点是能够利用纸笔得到简单有效而且准确解；缺点是能够求解方程数量太少。Matlab和Maple提供了求方程解析解函数，能够说对数学演算提供了不少方便。/10/106第6页2，方程求解之图形放大法图形最大有点就是直观，试想假如我们有了函数准确图形，那么曲线和x轴交点就是我们要求方程解。所以我们能够利用图形工具得到方程解。当然，计算机上图形不可能等同于函数真实图形，因为计算机上图形是曲线上部分点轨迹而不是全部，所以经过图形不可能得到方程准确解，甚至它只是一个比较粗略解，当然，经过对图形放大能够得到更准确一些解。同时，这种方法也不适应大量数据处理。/10/107第7页 方程 f(x)=0 1）建立坐标系，画曲线f(x)；2）观察曲线f(x)与x轴相交交点；3）将其中一个交点进行局部放大；4）该交点横坐标值就是方程根。2图形放大法/10/108第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，22内，所以将区间内，所以将区间-3-3，33缩缩小至小至-2-2，22，再观察！，再观察！该方程有几个根？欲寻找其中一个实根，而且到达一定精度。2图形放大法/10/109第9页 逐次缩小区间，观察一个根在-1.55-1.5之间。2图形放大法/10/1010第10页3,方程求解之迭代法迭代法理论以及方法出现，对方程求解有着里程碑式意义。其基本思想以下:需要求解方程：f(x)=0 （1）经过某种变形得：x=j(x)（2）从而求解方程（1）转化成为求解（2）得不动点。（满足条件x*=j(x*)点x*称为不动点）为得到方程不动点，能够结构迭代过程以下：xn+1=j(xn)，n=0，1，x0 定义为迭代初值。/10/1011第11页解：解：第一步第一步 结构迭代函数：x=j(x)例：例：用迭代方法求解方程 x3-x2-x-1=0。/10/1012第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/3);%j2(y)z=1+1/z+1/z2;%j3(z)endX，y，z/10/1013第13页序号j3(x)序号j2(x)j3(x)1j2(x)1.44223.000081.81751.813621.65371.444491.83851.855431.75322.1716101.83891.829441.79951.6725111.83911.845451.82091.9554121.83921.835561.83081.7730131.83921.841671.83541.8822j1(x)迭代是失败（迭代不收敛）。准确解：x=1.8393计算结果计算结果/10/1014第14页 迭代函数j2(x)和j3(x)选取是成功。准确解为 x=1.8393。而且选取函数j2(x)、j3(x)其收敛速度不一致，前者速度快些！结论1、当碰到迭代不收敛时有什么处理方法？2、怎样提升收敛速度？对于给定方程 f(x)=0,有各种方式将它改写成等价形式 x=j(x)。但主要是怎样改写使得序列收敛？怎样改写使得序列收敛？返 回/10/1015第15页当今最流行迭代法是牛顿法以及由此改进一些方法，比如拟牛顿法等。其基本思想就是结构迭代格式是利用函数导数，这类方法有收敛速度快，稳定性好等特点。对低维和高维情况都适合，也是当今一些软件均采取方法。当然，因为需要函数导数信息，所以自然对不可微问题受到制约。该方法迭代格式为：xk+1=xk-f(xk)/f(xk)/10/1016第16页例:利用牛顿法求方程x3-x2-x-1=0根.方法:第一步,给出函数导函数3x2-2x-1;第二步,给出函数迭代格式:xk+1=xk-f(xk)/f(xk);设置一定精度要求,到达即终止.定义函数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);end/10/1017第17页三，解方程函数格式及例子Matlab对方程求解提供了以下一些函数：（1）多项式求根；（2）线性方程组求解；（3）普通非线性方程（组）求解：/10/1018第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.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)roots()语句使用方法/10/1019第19页4、线性方程组：AX=bAX=b 其中A是mn阶矩阵，b是m维向量。x=A b or x=inv(A)*b特点：只能求出一个特解。/10/1020第20页solution to the following linear system of equations:You can formulate and solve the problem asA=3 11-2;1 1-2;1-1 1;b=7;4;19;x=Abx=13.2188 -2.3438 3.4375/10/1021第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()语句使用方法返 回/10/1022第22页函数fzero格式：fzero(函数名,初值或区间)Example 1.Calculate by finding the zero 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/10/1023第23页1、方程(组)，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软件直接求解法/10/1024第24页输出:1/2/a*(-b+(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）符号解/10/1025第25页例2：解方程：x3-2x2=x-1解：s=solve(x3-2*x2=x-1)double(s)2）数字解该方程是否有实根？vpa(s,10)solve()语句使用方法/10/1026第26页 例3 求解方程：tan(x)-sin(x)=0 3）无穷解输入:solve(tan(x)-sin(x)=0)输出：0 (不能给出全部解）(tx1.m)solve()语句使用方法/10/1027第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或 b=2solve()语句使用方法/10/1028第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)x=0.59905375664056731520568183824539y=2.3959314023778168490940003756591z=2.0050148409816158357003177860955输出：fsolve()语句使用方法/10/1029第29页解：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）运行结果：Optimization Terminated Successfully y=0.5990 2.3959 2.0050fsolve()语句使用方法/10/1030第30页fsolve()函数第三个输入是函数第三个输入是options,它是一个结构型数据它是一个结构型数据,能够经能够经过函数过函数optimset()进行设定进行设定.当不进行设定时采取缺省设置当不进行设定时采取缺省设置.fsolve函数还能够有后面参数设定函数还能够有后面参数设定,这一功效在有时候非常有用这一功效在有时候非常有用.fsolve()语句使用方法比如,求解函数sin(ax)-x=0最小正解.分析:已知知识我们知道,该方程没有解析形式解,也就是说极难得到这个解和a之间详细关系.数值方法能够对a先取确定值,这么变成一个一元方程,轻易进行求解,当我们改变a时,就得到很多这么解,经过这么方法,我们能够得到解和a之间一些大致关系./10/1031第31页详细方法以下详细方法以下:建立函数建立函数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(B)/10/1032第32页/10/1033例：某物体边缘呈圆形，经过测量边缘上是一个点坐标，数据为（见程序）。使用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.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第33页区间方法对于一个闭区间上连续函数，我们有一个0点存在定理，利用这个定理，能够对不可微函数求得函数0点。基本思想是经过判断函数在端点处函数值异号能够确定函数在开区间上最少有一个0点，然后经过缩小区间得到解近似。优点是不需要函数导数信息，而且只要有0点就一定能够得到；缺点是相对牛顿法等速度较慢。/10/1034第34页 问题关键在于圆桶到底能承受多大速度碰问题关键在于圆桶到底能承受多大速度碰撞撞?圆桶和海底碰撞时速度有多大圆桶和海底碰撞时速度有多大?工程师们进行了大量破坏性试验工程师们进行了大量破坏性试验,发觉圆发觉圆桶在直线速度为桶在直线速度为40 ft/s 冲撞下会发生破裂冲撞下会发生破裂,剩下剩下问题就是计算圆桶沉入问题就是计算圆桶沉入300 ft 深海底时深海底时,其末速其末速度终究有多大度终究有多大?问题分析问题分析引例分析和求解/10/1035第35页 1.使用使用55加仑圆桶加仑圆桶;(1加仑加仑=3.7854升升)2.装满放射性废物时圆桶重量为装满放射性废物时圆桶重量为 W=527.436磅磅 (1 磅磅=0.4526千克千克)3.在海水中圆桶受到浮力在海水中圆桶受到浮力 B=470.327磅磅 4.圆桶下沉时受到海水阻力圆桶下沉时受到海水阻力 D=C v C 为常数为常数,经测算得经测算得:C=0.08.5.建立坐标系建立坐标系,取垂直向下为坐标方向取垂直向下为坐标方向 y,海平面为坐标原点海平面为坐标原点.y0 问题假设问题假设引例/10/1036第36页依据牛顿第二定律依据牛顿第二定律,圆桶下沉时应满足微分方程圆桶下沉时应满足微分方程:建立模型建立模型引例/10/1037第37页 为了求出圆桶与海底碰撞速度为了求出圆桶与海底碰撞速度v(t),需要求出圆需要求出圆桶下沉到海底桶下沉到海底300英尺时时间英尺时时间 t,再计算再计算v(t)，要做到要做到这一点是十分困难这一点是十分困难.若将速度若将速度v 看成是海水深度看成是海水深度y 函函数数,即即由复合函数求导法知由复合函数求导法知 建立模型建立模型引例/10/1038第38页借助数值方法求出借助数值方法求出 v(300)=45.1ft/s，显然大于，显然大于40ft/s。结论：放射性废物不能随意放入公海！建立模型建立模型引例/10/1039第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、求解线性方程组：实 验 内 容/10/1040第40页3、求以下方程组解使用命令 solve()或 fsolve()。x0=-5,-5实 验 内 容/10/1041第41页4.炮弹发射视为斜抛运动，已知初始速度为 200 m/s，问要击中水平距离360m、垂直距离160m 目标，当忽略空气阻力时，发射角应多大？深入思索深入思索：假如要考虑水平方向阻力，且设阻力与（水平方向）速度成正比，系数为 0.1（1/s），结果又怎样？炮弹质量假设为10千克？实 验 内 容/10/1042第42页5.讨论:有这么电视节目,对一件商品要求观众猜价格,主持人对观众所猜数目标回答是高或者低(相对于实际价格).现在问题时,对于价格在1000元之内商品,寻找一个最好方法,确保在最少次数猜中商品价格.试给出竞猜5次得到正确答案概率位多大.(设商品最小价格单位为元)实 验 内 容结结结结 束束束束提醒:利用二分法思想,能够比较快得到竞猜次数最大值.一个更有趣方法是二叉树,其深度就是对应解.对于平均次数,能够采取模拟方法,比如产生10000次随机数,看每次猜中次数,平均就得到平均次数./10/1043第43页