实验5非线性方程求根及其MATLAB实现

分类数学实验之方程求近似实根问题、高射炮的控制区域问题对应书本的第二章数值计算问题的2.1及2.4考虑方程f(x)=0⑴求根分为两步：（１）先确定某个根的近似值；（２）再将初始近似值加工成满足精度要求的结果。１.两分法迭代（理论基础：零点定理）设f(x)∈C[a,b]，f(a)·f(b)0。区间(a,b)就是方程⑴根的存在区间，再用下面的方法改善根的精度。方程求根数值算法的基本思想取[a,b]的中点x0=(a+b)/2，若f(x0)=0，则x0即是根；否则，f(a)·f(x0)0，令a1=a，b1=x0(取[a,b]的左半部)；f(x0)·f(b)0，令a1=x0，b1=b(取[a,b]的右半部)。aba+b2f(a)f((a+b)/2)f(b)f(x)xy将上述做法重复n次，得到n个小区间，且bn-an=(b-a)/2n，[a,b]经过一次这样对原区间[a,b]的处理，得到了一个新的有根区间[a1,b1]，且b1-a1=(b-a)/2，[a,b][a1,b1]。１.两分法迭代当n足够大时即可达到满意的精度σ。[a1,b1]…[an,bn]。需要的迭代次数为：ln(b-a)-lnN=ln2例1求方程x3+1.1x2+0.9x–1.4=0的一个实根。使误差不超过10-3１.两分法迭代解：（１）首先观测图形，作f(x)的图像:ezplot('x^3+1.1*x^2+0.9*x-1.4')x=-1:0.01:1;y=x.^3+1.1.*x.^2+0.9.*x-1.4;%函数表达式figure;plot(x,y,'LineWidth',2)%画出图形holdon;y1=zeros(size(x));%y1＝0plot(x,y1,'r','LineWidth',4);先绘制其图形观察一下：bisect.m（２）按两分法的思想，进行迭代求根。为了使得程序具有通用性，将方程的表达式写成一个函数：functiony=myequation(x);y=x.^3+1.1.*x.^2+0.9.*x-1.4;这样修改y的表达式，即可求出其他方程的实根。（3）运行二分法的程序：(该程序在matlab中打开)x0=-1;x1=2;%初始值根始终在x0和x1之间i=0;delt=10e-5;%精度N=ceil((log(x1-x0)-log(delt))/log(2));%理论迭代次数whileabs(x1-x0)delt%达不到我们需要的精度y0=myequation(x0);y1=myequation(x1);%求两个初始点处的y值x2=(x0+x1)./2;%二分y2=myequation(x2);%求在中点处的值ify2==0%恰好是方程的根x0=x2;x1=x2;elseify0*y20%方程的根落在x0和x2之间x1=x2;%新的两个端点为：x0=x0;x1=x2。else%方程的根落在y1和y2之间x0=x2;%新的两个端点为：x0=x2;x1=x1。endi=i+1;enddisp(‘方程的根为：x=’);x0disp(‘实际迭代次数为：i=');idisp('理论迭代次数为：N=');Nbisect.m两分法迭代的加速：aba+b2f(a)f((a+b)/2)f(b)f(x)xyc以f(a),f(b)的连线在x轴的交点c作为新的出发点。程序：fastbisect.m()-()c()()afbbfacfbfa点坐标可求出为：注：此方法不一定能加速!!!减速举例：aba+b2f(a)f(b)f(x)xyc程序：fastbisect.mf((a+b)/2)加速情况举例：２.不动点迭代称满足方程f(x)=x的点x为函数f的不动点.求函数f的不动点。可以从一个初始点x0出发，以格式xn+1=f(xn)进行迭代；x1=f(x0)，x2=f(x1)，…，xk+1=f(xk)，…得到x0,x1,x2,……,xn,…..如果该数列是收敛的，则11*)lim(*lim;()nnnnxfxxxf将方程f(x)=0⑴化为等价方程x=φ(x)（2）取某个定数x0，做数列{xn}，其中，x1=φ(x0)，x2=φ(x1)，…，xk+1=φ(xk)，…（3）设φ(x)连续，且axnnlima就是方程f(x)=0的根。a就是函数φ的一个不动点，即a=φ(a)等价于f(a)=0求根原理：y=f(x)x0f(x0)x1=f(x0)f(x1)x2=f(x1)y=xx收敛的迭代：观察动态演示：发散的迭代：观察动态演示：x0kx)(xy由以上图形可以看出，序列的收敛速度，取决于曲线在根附近的斜率)(')()(1111kkkkkkkkkxxxxxxxx所以在根x*附近，’(x)恒小于1，则此迭代序列收敛，若’(x)1,则此序列发散。２.不动点迭代例2求方程0142xx在x=3附近的近似实根。解：可将方程写成下三种形式：x=14–x2,114xx12142xxxxxfunctionf=iterfun(x)%f=14-x.^2;%f=14./(x+1);f=x-(x.^2+x-14)./(2*x+1);(1)将三种迭代形式写成函数存起来：(2)演示程序使用的程序：iter.m２.不动点迭代收敛性蛛网图一般地，若函数φ(x)在含不动点ξ的某邻域内一阶导数连续，且1|)(|00.20.40.60.8100.10.20.30.40.50.60.70.80.91My=x00.20.40.60.8100.10.20.30.40.50.60.70.80.91My=x则存在一个邻域Δ：|ξ-x|δ，对任何的x0∈Δ，其迭代序列必收敛。同时我们还提供一个高级程序：zhuwang.m，绘制出迭代格式为41xx的蛛网图。此处用到一个函数:quiver,画出矢量方向。quiver(x,y,u,v):画出在(x,y)处，方向为u,v的矢量供学有余力的、感兴趣的同学们好好体会学习，体会Matlab绘图的精妙。3.牛顿迭代法记[a,b]为方程f(x)=0的根的存在区间，f(a)与f(b)异号，且对于每个x∈[a,b]，f'(x)≠0,f(x)保持符号不变。取x0∈[a,b]，对f(x)用微分中值定理，近似地，有f(x)≈p(x)=f(x0)+f'(x0)(x-x0)令p(x)=0，得到f(x)=0的近似根000()()fxxxfx记为x1=为改善根的精确程度，反复实施这一过程，得到牛顿迭代公式:1()()kkkkfxxxfxxk+1就是从线性方程f(xk)+f'(xk)(x–xk)=0(1)中解得的根x.以点(xk,f(xk))为切点，曲线y=f(x)的切线方程恰为：y=f(xk)+f'(xk)(x–xk)(2)方程（2）与y=0相联立，得到与x轴的交点，该交点即xk+1x*x0x0,f(x0)x10100()()fxxxfxx1,f(x1)x2例3求方程3310xx在x0=2附近的近似实根。准确到小数点后4位数字解：3()31fxxx2'()33fxx迭代公式为：3123133kkkkkxxxxx计算步骤如下：（1）选x0=2,按照迭代公式计算x1；（2）若|x1-x0|=0.00001,终止迭代；否则，x0=x1;转（1）；（3）输出迭代次数和近似根编程步骤如下：（1）写函数：function[fun,dfun]=fun0(x)fun=x^3-3*x-1;%求原函数的值dfun=3*x^2-3;%求一阶导数的值；（2）写牛顿迭代法的程序：x0=2;[fun,dfun]=fun0(x0);x1=x0-fun/dfun;i=1;whileabs(x1-x0)0.00001x0=x1;[fun,dfun]=fun0(x0);x1=x0-fun/dfun;i=i+1;enddisp('thesolutionisx1=')x1disp('theitertimeis')inewtuniter.m第二部分：高射炮的控制区域问题问题描述：高射炮发射的炮弹在空中呼啸而过划出一条抛射线。设坐标原点为发射点，x(t)和y(t)分别是高射炮在t时刻的空间坐标。用v0为炮弹出膛时的初速度，α为高射炮发射角度，g是重力加速度。问：（1）试推高射炮抛射的弹道曲线的参数方程；（2）当炮弹出膛速度v0确定时，它的最远射程？（3）当炮击目标确定后，如何调整发射角度，使炮弹能准确地落在目标位置处爆炸？（4）一门高射炮可以控制什么样的空间区域？第二部分：高射炮的控制区域问题问题描述：高射炮发射的炮弹在空中呼啸而过划出一条抛射线，设坐标原点为发射点，抛射的弹道曲线的参数方程为：问题：（1）当炮弹出膛速度v0确定时，它的最远射程？（2）当炮击目标确定后，如何调整发射角度，使炮弹能准确地落在目标位置处爆炸？（3）一门高射炮可以控制什么样的空间区域？其中v0为炮弹出膛时的初速度，α为高射炮发射角度，g是重力加速度。初值条件：00(0)0,cos(0)0,sindxxdtdyygtdt即可导出弹道曲线的参数方程为(1)模型建立与分析消去参数t，可以得到弹道曲线的直角坐标的表达式：2220sec(tan)2gxyx这是我们十分熟悉的抛物线。（1）（2）就形成一族抛物线段.所以，（1）或（2）可以看作是带有参数α的平面曲线族。当发射角变化时，取900，相应的弹道曲线(2)控制区与包络线包络曲线：如果一条曲线的每一点都与曲线簇中的某一曲线相切，则称这条曲线为该曲线簇的包络。对曲线簇F(x,y,c)=0中，对参数c求导，与它自己联立得到方程组：由这个方程组确定的曲线x=x(c)即y=y(c)称为单参数曲线簇F(x,y,c)的包络。包络曲线与x=0及y=0所形成的区域，即高射炮的控制区域。0),,(0),,('cyxFcyxFc(额外知识：)如果曲线族的参数方程形式已给出，则类似可推出其包络线应满足的判别方程组为：（由于此部分超出同学们知识范围，过程略）(,)(,)0xxtyytxyxytt(额外知识：)2202022gyxg可得包络曲线的表达式：3.模型求解：首先建立平面直角坐标系，设坐标原点为发射点。初始速度向量为：发射点到落点的距离是射程，初始速度不变时，落点仅与发射角有关。落点时，y=0,于是有：y=解之，得落点所对应的参数值：02sinvtg此时，落地点20002sincossin2vvxvgg显然，若要使射程最大，即x最大，则发射角度4模型验证：取角度间隔,10060v为方便计，固定初始速度对于不同的发射角度，分别计算出它们的射程和空中轨迹，并画出射程最大时的空中轨迹。请查看matlab源程序：paodan.m模型的拓展：（1）若将此问题拓展到铅球投掷试验中，应做什么修改？（2）若已知目标位置，如何调整发射角？