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

bisect.m
两分法迭代的加速：
y 以f(a),f(b)的连线在x轴的交点c作为新的出
af (b) - bf (a) c点坐标可求出为：c 发点。 f (b) f (a)
注：此方法不一定能加速!!! f(x) f(b) 减速举例： c
a
a+b 2
b
x
f(a)
f((a+b)/2)
程序：fastbisect.m
f ( x0 ) 记为x1= x x0 f ( x 0 )
为改善根的精确程度，反复实施这一过程，得到牛顿 迭代公式: f (x )
xk 1 xk f ( xk )
k
xk+1就是从线性方程 f (xk) + f '(xk)(x – xk) = 0 中解得的根x.
(1)
以点(xk , f (xk))为切点，曲线y = f (x)的切线方程恰为： y = f (xk) + f '(xk)(x – xk) (2) 方程（2）与 y=0相联立，得到与x轴的交点，该交点 即xk+1
x x(t , ) y y (t , ) x y x y t t 0
g 2 可得包络曲线的表达式： y 2x 2 g 20
02
3. 模型求解： 首先建立平面直角坐标系，设坐标原点为发射点。 初始速度向量为： 发射点到落点的距离是射程，初始速度不变时， 落点仅与发射角有关。落点时，y=0,于是有： y= 解之，得落点所对应的参数值：
xk 3 3 xk 1 迭代公式为：xk 1 xk 3 xk 2 3
计算步骤如下：
（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;%求原函数的值 （2）写牛顿迭代法的程序： dfun=3*x^2-3;% 求一阶导数的值； newtuniter.m x0=2;
第二部分：高射炮的控制区域问题
问题描述：
高射炮发射的炮弹在空中呼啸而过划出一条抛射线， 设坐标原点为发射点，抛射的弹道曲线的参数方程为：
其中v0为炮弹出膛时的初速度，α为高射炮发射角度， g是重力加速度。
问题： （1）当炮弹出膛速度v0确定时，它的最远射程？ （2）当炮击目标确定后，如何调整发射角度， 使炮弹能准确地落在目标位置处爆炸？
y
f(b)
f(x)
a
a+b 2
b x
f(a) f((a+b)/2)
１.两分法迭代
经过一次这样对原区间[a, b]的处理，得到了一个新 的有根区间[a1 ,b1]，且 b1-a1 = (b-a)/2， [a, b] [a1 ,b1]。 将上述做法重复n次，得到n个小区间，且 bn-an = (b-a)/2n ， … [an ,bn]。 [a, b] [a1 ,b1]
F ( x, y , c ) 0 ' Fc ( x, y , c) 0
由这个方程组确定的曲线x=x(c) 即y=y (c) 称为单参数 曲线簇F(x,y,c)的包络。包络曲线与x=0及y=0所形成的 区域，即高射炮的控制区域。
(额外知识：) 如果曲线族的参数方程形式已给出，则类似可推出 其包络线应满足的判别方程组为：（ 由于此部分超出 同学们知识范围，过程略）
x=-1:0.01:1; y=x.^3+1.1.*x.^2+0.9.*x-1.4; %函数表达式 figure; plot(x,y,'LineWidth',2) %画出图形 hold on; y1=zeros(size(x)); %y1＝0 plot(x,y1,'r','LineWidth',4);
[fun,dfun]=fun0(x0); x1=x0-fun/dfun;i=1; while abs(x1-x0)>0.00001 x0=x1; [fun,dfun]=fun0(x0); x1=x0-fun/dfun; i=i+1; end disp('the solution is x1=') x1 disp('the iter time is ') i
供学有余力的、感兴趣的同学们好好体会学习，体会 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的近似根
x0=-1; x1=2;%初始值 根始终在x0和x1之间 i=0; delt=10e-5;%精度 N=ceil((log(x1-x0)-log(delt))/log(2));%理论迭代次数 while abs(x1-x0)>delt %达不到我们需要的精度 y0=myequation(x0); y1=myequation(x1);%求两个初始点处的y值 x2=(x0+x1)./2; %二分 y2=myequation(x2); %求在中点处的值 if y2==0 %恰好是方程的根 x0=x2; x1=x2; elseif y0*y2<0 % 方程的根落在x0和x2之间 x1=x2; %新的两个端点为：x0=x0;x1=x2。 else % 方程的根落在y1和y2之间 x0=x2; %新的两个端点为：x0=x2;x1=x1。 end i=i+1; end disp(‘方程的根为：x=’);x0 disp(‘实际迭代次数为：i=');i disp('理论迭代次数为：N=');N
(2)演示程序 使用的程序：iter.m
２.不动点迭代
收敛性
蛛网图
一般地，若函数φ(x)在含不动点ξ的某邻域内一阶导 数连续，且 | ( ) | 1 则存在一个邻域Δ：|ξ-x|<δ，对任何的x0∈Δ， 其迭代序列必收敛。
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2
x0
２.不动点迭代
例2 求方程
x x 14 0
2
在 x = 3附近的近似实 根。
解：可将方程写成下三种形式： x = 14 – x 2 ,
14 x x 1
x x 14 x x 2x 1
2
(1)将三种迭代形式写成函数存起来：
function f=iterfun(x) % f=14-x.^2; % f=14./(x+1); f=x-(x.^2+x-14)./(2*x+1);
求根 原理：
（ 3）
lim x n a
n
a就是函数φ的一个不动点，即 a =φ(a) 等价于 f (a ) = 0
a就是方程f (x) = 0的根。
f(x0)
收敛的迭代：观察动态演示： y=x
f(x1)
x0 x2= f(x1) x
y=f(x)
x1= f(x0)
发散的迭代：观察动态演示：
gx 2 sec2 y (tan ) x 2 20
（2）
这是我们十分熟悉的抛物线。
(2)控制区与包络线
当发射角变化时，取 0 90 ，相应的弹道曲线 就形成一族抛物线段.所以，（1）或（2）可以看作是 带有参数α的平面曲线族。
(额外知识：)
包络曲线：如果一条曲线的每一点都与曲线簇中的某 一曲线相切，则称这条曲线为该曲线簇的包络。 对曲线簇F(x,y,c)=0中，对参数c求导，与它自己联立 得到方程组：
（２）按两分法的思想，进行迭代求根。
为了使得程序具有通用性，将方程的表达式写成 一个函数：
function y=myequation(x);
y=x.^3+1.1.*x.^2+0.9.*x-1.4;
这样修改y的表达式，即可求出其他方程的实根。
（3）运行二分法的程序： bisect.m (该程序在matlab中打开)
１.两分法迭代（理论基础：零点定理）
取[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]的右半部)。
分类数学实验之
方程求近似实根问题、 高射炮的控制区域问题
对应书本的第二章 数值计算问题的2.1 及2.4
方程求根数值算法的基本思想 考虑方程 f (x) = 0 ⑴ 求根分为两步： （１）先确定某个根的近似值； （２）再将初始近似值加工成满足精度要求的结果。 设f (x)∈C[a, b]，f (a)· f (b) <0。区间(a, b)就是方程⑴ 根的存在区间，再用下面的方法改善根的精度。
第二部分：高射炮的控制区域问题
问题描述： 高射炮发射的炮弹在空中呼啸而过划出一条抛射线。 设坐标原点为发射点，x(t)和y(t)分别是高射炮在t时 刻的空间坐标。用v0为炮弹出膛时的初速度，α为高 射炮发射角度，g是重力加速度。问： （1）试推高射炮抛射的弹道曲线的参数方程； （2）当炮弹出膛速度v0确定时，它的最远射程？ （3）当炮击目标确定后，如何调整发射角度， 使炮弹能准确地落在目标位置处爆炸？ （4）一门高射炮可以控制什么样的空间区域？
x1 x0
f ( x0 ) f ( x0 )
x0,f(x0)
x1,f(x1)
x*
x2
x1
x0
例3 求方程
x3 3 x 1 0
在 x0 = 2附近的近似实根。 准确到小数点后4位数字 解： f ( x ) x 3 3 x 1
f '( x ) 3 x 2 3
加速情况举例：
y
f((a+b)/2)
f(b)
f(x)
a
c
a+b 2
b
x
f(a)
程序：fastbisect.m
２.不动点迭代

称满足方程 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,….. 如果该数列是收敛的，则
0.1 1
y=x
0.9 0.8
y=x M
M
0.7 0.6 0.5 0.4 0.3 0.2
0.1
0
0来自百度文库
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
同时我们还提供一个高级程序：zhuwang.m，绘制出 迭代格式为
4 x 1 x
的蛛网图。
此处用到一个函数:quiver,画出矢量方向。
quiver(x,y,u,v):画出在(x,y)处，方向为u,v的矢量
由以上图形可以看出，序列
xk 的收敛速度，取决于曲线 y ( x) 在根附近的斜率
xk 1 xk ( xk ) ( xk 1 ) ' ( k ) xk xk 1 xk xk 1
所以在根x*附近，’(x) 恒小于1 则此迭代序列收敛， 若’(x) 1,则此序列发散。
（3）一门高射炮可以控制什么样的空间区域？
(1)模型建立与分析
x(0) 0, 初值条件： y (0) 0, dx 0 cos dt dy 0 sin gt dt
即可导出弹道曲线的参数方程为 （1）
消去参数t，可以得到弹道曲线的直角坐标的表达式：
x* lim xn 1;
n
f ( x*) lim f ( xn 1 )
n
将方程 f (x) = 0 ⑴ 化为等价方程 x =φ(x) （ 2） 取某个定数x0，做数列{xn}，其中， x1 =φ(x0)，x2 =φ(x1)，…，xk+1 =φ(xk)，… 设φ(x)连续，且


当n足够大时即可达到满意的精度σ。 需要的迭代次数为：N= ln(b-a)-ln ln2
１.两分法迭代
先绘制其图形观察一下：
ezplot('x^3 + 1.1*x^2 + 0.9*x- 1.4') 例1 求方程 x 3 + 1.1x 2 + 0.9x – 1.4 = 0 的一个实根。使误差不超过10-3 解：（１）首先观测图形，作f(x)的图像: