数学实验第三次讲稿

数学实验第三次讲稿第1页/共44页数学实验第三讲－－方程求解2023/4/112第2页/共44页

1)掌握方程求解的三种解法：解析法、数值解法以及图形表示解的方法;2)学会使用MATLAB软件求解析解、数值解和图形解；

3）通过范例学习怎样建立方程模型和分析问题的思想。实验目的2023/4/113第3页/共44页一，建立方程例子【问题背景】一段时间,美国原子能委员会是按以下方式处理浓缩放射性废物的.他们将废物装入密封性能很好的圆桶中,然后扔到水深300英尺的海里.这种做法是否会造成放射性污染,很自然地引起了生态学家及社会各界的关注.原子能委员会一再保证,圆桶非常坚固,决不会破漏,这种做法是绝对安全的.然而一些工程师们却对此表示怀疑,他们认为圆桶在海底相撞时有可能发生破裂.由此双方展开了一场笔墨官司.究竟谁的意见正确呢?只能让事实说话了!2023/4/114第4页/共44页数学建模的一般步骤模型准备模型假设模型构成模型检验模型分析模型求解模型应用2023/4/115第5页/共44页二，方程求解1，解析方法2，图形放大法3，迭代方法4，区间方法2023/4/116第6页/共44页1，方程求解之解析方法主要针对一些比较简单的方程以及方程组，比如多项式方程等。同学们以前对方程的求解也是针对这样一些方程进行的。该方法的优点是可以利用纸笔得到简单有效并且精确的解；缺点是可以求解的方程数量太少。Matlab和Maple提供了求方程解析解的函数，可以说对数学演算提供了不少方便。2023/4/117第7页/共44页2，方程求解之图形放大法图形最大的有点就是直观，试想如果我们有了函数准确的图形，那么曲线和x轴的交点就是我们要求的方程的解。因此我们可以利用图形工具得到方程的解。当然，计算机上的图形不可能等同于函数的真实图形，因为计算机上的图形是曲线上部分点的轨迹而不是全部，因此通过图形不可能得到方程的精确解，甚至它只是一个比较粗略的解，当然，通过对图形的放大可以得到更精确一些的解。同时，这种方法也不适应大量的数据处理。2023/4/118第8页/共44页

方程

f(x)=01）建立坐标系，画曲线f(x)；

2）观察曲线f(x)与x轴相交的交点；

3）将其中一个交点进行局部放大；

4）该交点的横坐标值就是方程的根。2．图形放大法2023/4/119第9页/共44页例:

求方程

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/4/1110第10页/共44页

逐次缩小区间，观察一个根在-1.55~-1.5之间。2．图形放大法2023/4/1111第11页/共44页3,方程求解之迭代法迭代法的理论以及方法的出现，对方程求解有着里程碑式的意义。其基本思想如下:

需要求解的方程：

f(x)=0（1）经过某种变形得：x=j(x)

（2）从而求解方程（1）转化成为求解（2）得不动点。（满足条件x*=j(x*)的点x*称为不动点）为得到方程的不动点，可以构造迭代过程如下：

xn+1=j(xn)，n=0，1，…x0定义为迭代初值。2023/4/1112第12页/共44页解：

第一步构造迭代函数：

x=j(x) 例：用迭代方法求解方程

x3

-x2-x-1=0。2023/4/1113第13页/共44页第二步迭代

设定初值x0=1，

xn+1=j(xn)，n=0，1，…

用MATLAB编程（died2.m文件)x=1;y=1;z=1;(初始点）fork=1:20x=x^3-x^2-1;%j1(x)y=(y^2+y+1)^(1/3);%j2(y)z=1+1/z+1/z^2;%j3(z)endX，y，z2023/4/1114第14页/共44页序号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计算结果2023/4/1115第15页/共44页

迭代函数j2(x)和j3(x)的选取是成功的。精确解为

x=1.8393。并且选取函数j2(x)、j3(x)其收敛速度不一致，前者的速度快些！结论1、当遇到迭代不收敛时有什么解决办法？2、如何提高收敛速度？

对于给定的方程

f(x)=0,有多种方式将它改写成等价的形式

x=j(x)。但重要的是如何改写使得序列收敛？返回2023/4/1116第16页/共44页当今最流行的迭代法是牛顿法以及由此改进的一些方法，比如拟牛顿法等。其基本的思想就是构造迭代格式是利用函数的导数，这类方法有收敛速度快，稳定性好等特点。对低维和高维情况都适合，也是当今一些软件均采用的方法。当然，因为需要函数的导数信息，所以自然对不可微的问题受到制约。该方法的迭代格式为：

xk+1=xk-f(xk)/f’(xk)2023/4/1117第17页/共44页例:利用牛顿法求方程x3-x2-x-1=0的根.方法:第一步,给出函数的导函数3x2-2x-1;

第二步,给出函数的迭代格式:

xk+1=xk-f(xk)/f’(xk);

设置一定的精度要求,达到即终止.定义函数m文件:functionff=mynewton(x)ff=(x^3–x^2-x-1)/(3*x^2-2*x-1);定义命令m文件:x0=1;x1=x0-mynewton(x0);whileabs(x1-x0)>0.0001x0=x1;x1=x0-mynewton(x0);end2023/4/1118第18页/共44页三，解方程函数格式及例子Matlab对方程的求解提供了以下的一些函数：（1）多项式求根；（2）线性方程组求解；（3）一般的非线性方程（组）求解：2023/4/1119第19页/共44页输出:

-1.2131-0.9017+0.5753i-0.9017-0.5753i-0.2694+0.9406i-0.2694-0.9406i0.4168+0.8419i0.4168-0.8419i0.8608+0.3344i0.8608-0.3344i例：求解多项式方程x9+x8+1=0输入:p=[1,1,0,0,0,0,0,0,0,1];roots(p)roots()语句的用法2023/4/1120第20页/共44页4、线性方程组：AX=b

其中A是m×n阶矩阵，b是m维向量。x=A\borx=inv(A)*b特点：只能求出一个特解。2023/4/1121第21页/共44页solutiontothefollowinglinearsystemofequations:YoucanformulateandsolvetheproblemasA=[311-2;11-2;1-11];b=[7;4;19];x=A\bx=13.2188-2.34383.43752023/4/1122第22页/共44页例:AX=b,

解:输入:A=[123;456;789];b=[6;14;-3];x1=A\b,x2=inv(A)*b输出: 警告:系统的秩不足.解不唯一.1、题中rank(A)=rank(A|b)=2<3,该方程组有无穷解。2、输出结果是否一致？3、如何求方程组的全部解？思考A\b和inv()语句的用法返回2023/4/1123第23页/共44页函数fzero格式：fzero(‘函数名’,初值或区间)Example1.Calculatebyfindingthezeroofthesinefunctionnear3.x=fzero(@sin,3)x=3.1416Example2.Tofindthezeroofcosinebetween1and2x=fzero(@cos,[12])x=1.57082023/4/1124第24页/共44页1、方程(组)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)2、方程(组)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)fun.mfunctionf=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软件直接求解法2023/4/1125第25页/共44页输出:[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]①单变量方程solve()语句的用法例1：求解方程ax2+bx+c=0输入: x=

solve('a*x^2+b*x+c')或solve('a*x^2+b*x+c=0')1）符号解2023/4/1126第26页/共44页例2：解方程：x3-2x2=x-1解：

s=solve('x^3-2*x^2=x-1')double(s)2）数字解该方程是否有实根？vpa(s,10)solve()语句的用法2023/4/1127第27页/共44页

例3求解方程：tan(x)-sin(x)=0

3）无穷解输入:solve('tan(x)-sin(x)=0')输出：

0(不能给出全部解）

(tx1.m)solve()语句的用法2023/4/1128第28页/共44页输入：[x,y]=solve('x^2*y^2,x-(y/2)-b')输出：x=[0]，y=[-2*b][0]，[-2*b]（符号解）

[b]，[0][b]，[0]v=[x,y]②多变量方程组例4或b=2solve()语句的用法2023/4/1129第29页/共44页例6：求解方程组解①输入：[x,y,z]=solve('sin(x)+y^2+log(z)-7=0','3*x+2^y-z^3+1=0','x+y+z-5=0','x','y','z')x=0.59905375664056731520568183824539y=2.3959314023778168490940003756591z=2.0050148409816158357003177860955输出：fsolve()语句的用法2023/4/1130第30页/共44页解②：1）建立方程组的M-函数文件(nxxf.m)

functioneq=nxxf(x)eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;eq(2)=3*x(1)+2^x(2)-x(3)^3+1;eq(3)=x(1)+x(2)+x(3)-5;运行程序(test4.m)y=fsolve('nxxf',[1,1,1])3）运行结果：OptimizationTerminatedSuccessfullyy=0.59902.39592.0050fsolve()语句的用法2023/4/1131第31页/共44页fsolve()函数的第三个输入是options,它是一个结构型数据,可以通过函数optimset()进行设定.当不进行设定时采用缺省设置.fsolve函数还可以有后面的参数设定,这一功能在有的时候非常有用.fsolve()语句的用法例如,求解函数sin(ax)-x=0的最小正解.分析:已知的知识我们知道,该方程没有解析形式的解,也就是说很难得到这个解和a之间的具体关系.数值的方法可以对a先取确定的值,这样变成一个一元方程,容易进行求解,当我们变化a时,就得到很多这样的解,通过这样的方法,我们可以得到解和a之间的一些大致关系.2023/4/1132第32页/共44页具体方法如下:建立函数m文件:functionff=funpara(x,a)ff=sin(a*x)-x;fsolve()语句的用法建立相应的命令m文件:B=zeros(100,1);fora=1:100x0=pi/(2*a)+0.01;B(a)=fsolve(@funpara,x0,[],a);endplot(B)2023/4/1133第33页/共44页2023/4/1134例：某物体的边缘呈圆形，通过测量边沿上是一个点的坐标，数据为（见程序）。使用fsolve计算物体边缘的曲线方程。解：假定物体边缘的曲线方程为（x-a）^2+(y-b)^2=r^2为了得到a，b和r的值。将上述十一组数据代人方程，得到（xi-a）^2+(yi-b)^2=r^2（i=1,…,11）。次数方程个数是11个，未知量个数是3个。使用fsolve求解，其函数m文件为：functiony=funcircle(x)A=[6.76305.13132.4713-0.3435-2.3887-2.9927-1.95720.37783.24555.70426.9465];B=[23.287925.649226.726826.166824.153121.347018.669917.001016.888318.368820.9564];fori=1:11y(i)=(A(i)-x(1))^2+(B(i)-x(2))^2-x(3)^2;end第34页/共44页区间方法对于一个闭区间上的连续函数，我们有一个0点存在定理，利用这个定理，可以对不可微的函数求得函数的0点。基本思想是通过判断函数在端点处的函数值异号可以确定函数在开区间上至少有一个0点，然后通过缩小区间得到解的近似。优点是不需要函数的导数信息，而且只要有0点就一定可以得到；缺点是相对牛顿法等速度较慢。2023/4/1135第35页/共44页

问题的关键在于圆桶到底能承受多大速度的碰撞?圆桶和海底碰撞时的速度有多大?

工程师们进行了大量破坏性的实验,发现圆桶在直线速度为40ft/s的冲撞下会发生破裂,剩下的问题就是计算圆桶沉入300ft深的海底时,其末速度究竟有多大?问题分析引例的分析和求解2023/4/1136第36页/共44页1.使用55加仑的圆桶;(1加仑=3.7854升)2.装满放射性废物时的圆桶重量为

W=527.436磅(1磅=0.4526公斤)3.在海水中圆桶受到的浮力B=470.327磅

4.圆桶下沉时受到海水的阻力D=Cv