数值分析实验报告97389

数值分析实验报告

（第二章）

实验题目：

分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程

f(x)=(x2+1)(x?1)5=0

的根x=1，观察不同初始值下的收敛性，并给出结论。

问题分析：

题目有以下几点要求：

1.不同的迭代法计算根，并比较收敛性。

2.选定不同的初始值，比较收敛性。

实验原理：

各个迭代法简述

二分法：取有根区间[a,b]的重点x0，确定新的有根区间[a1,b1]的区间长度仅为[a,b]区间长度的一版。对压缩了的有根区间[a1,b1]重复以上过程，又得到新的有根区间[a2,b2]，其区间长度为[a1,b1]的一半，如此反复，……，可得一系列有根区间，区间收敛到一个点即为根。

牛顿迭代法：不动点迭代法的一种特例，具有局部二次收敛的特性。迭代格式为

x n+1=x n?f(x n)

f′(x n)

,n=0,1,2,…

割线法：是牛顿法的改进，具有超线性收敛的特性，收敛阶为1.618. 迭代格式为

x n+1=x n?

f(x n)

(n)(n?1)

(x n?x n?1),n=1,2,…

史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为

y n=φ(x n)

z n=φ(y n)

x n+1=x n?

(y n?x n)2 z n?2y n+x n

这里φ(x)可采用牛顿迭代法的迭代函数。实验内容：

1.写出该问题的f(x)函数代码如下：

function py= f(x) syms k;

y=(k^2+1)*(k-1)^5;

yy=diff(y,k);

py(1)=subs(y,k,x);

py(2)=subs(yy,k,x); end

2.分别写出各个迭代法的迭代函数代码如下：

二分法：

function y=dichotomie(a,b,e) i=2;

m(1)=a;

while abs(a-b)>e

t=(a+b)/2;

s1=f(a);

s2=f(b);

s3=f(t);

if s1(1)*s3(1)<=0

b=t;

else

a=t;

end

m(i)=t;

i=i+1;

end

y=[t,i+1,m];

end

牛顿迭代法：

function

y=NewtonIterative(x,e)

i=2;

en=2*e;m(1)=x;

while abs(en)>=e

s=f(x);

t=x-s(1)/s(2);

en=t-x;

x=t;

m(i)=t;

i=i+1;

end

y=[x,i+1,m];

end

牛顿割线法：

function y=Secant(x1,x2,e) i=3;

m(1)=x1,m(2)=x2;

while abs(x2-x1)>=e

s1=f(x1);

s2=f(x2);

t=x2-(x2-x1)*s2(1)/(s2(1)-s1( 1));

x1=x2;

x2=t;

m(i)=t;

i=i+1;

end

y=[x2,i+1,m];

end

史蒂芬森迭代法：

Function p=StephensonIterative (x,e)

i=2;

m(2)=x;

en=2*e;

while abs(en)>=e

y=fai(x); z=fai(y);

t=x-(y-x)^2/(z-2*y+x); en=t-x;

x=t;

m(i)=t;

i=i+1;

end

p=[x,i+1,m];

end

3.因为φ(x)经常被使用，故可以写一个φ(x)函数。

代码如下：

function y=fai(x)

s=f(x);

y=x-s(1)/s(2);

end

4.可以绘制不同的图形来比较不同迭代法的收敛性和不同初值下的收敛性。代码如下：

clear all;

%相同初始值，不同迭代法下的收敛

x1=dichotomie(0,3,1e-10);

x2=NewtonIterative(0,1e-10);

x3=Secant(0,2,1e-10);

x4=StephensonIterative(0,1e-10);

[x1(2),x2(2),x3(2),x4(2)]

figure,

subplot(2,2,1),plot(x1(3:x1(2))),title('二分法');

subplot(2,2,2),plot(x2(3:x2(2))),title('牛顿迭代法');

subplot(2,2,3),plot(x3(3:x3(2))),title('牛顿割线法');

subplot(2,2,4),plot(x4(3:x4(2))),title('史蒂芬森迭代法');

figure,

subplot(2,2,1),plot((x1(4:x1(2)-1)-x1(1))./(x1(3:x1(2)-2)-x1(1))),tit le('二分法');

subplot(2,2,2),plot((x2(4:x2(2)-1)-x2(1))./(x2(3:x2(2)-2)-x2(1))),tit le('牛顿迭代法');

subplot(2,2,3),plot((x3(4:x3(2)-1)-x3(1))./(x3(3:x3(2)-2)-x3(1))),tit

le('牛顿割线法');

subplot(2,2,4),plot((x4(4:x4(2)-1)-x4(1))./(x4(3:x4(2)-2)-x4(1))),tit le('史蒂芬森迭代法');

%不同初始值，相同迭代法下的收敛性

x5=dichotomie(-1,1,1e-10);

x6=dichotomie(-2,3,1e-10);

x7=dichotomie(0,4,1e-10);

x8=dichotomie(-4,4,1e-10);

x9=NewtonIterative(-2,1e-10);

x10=NewtonIterative(-4,1e-10);

x11=NewtonIterative(4,1e-10);

x12=NewtonIterative(6,1e-10);

figure,

subplot(1,2,1),

plot(1:x1(2)-2,x1(3:x1(2)),1:x5(2)-2,x5(3:x5(2)),1:x6(2)-2,x6(3:x6(2) ),1:x7(2)-2,x7(3:x7(2)),1:x8(2)-2,x8(3:x8(2))),title('二分法');

subplot(1,2,2),

plot(1:x2(2)-2,x2(3:x2(2)),1:x9(2)-2,x9(3:x9(2)),1:x10(2)-2,x10(3:x10 (2)),1:x11(2)-2,x11(3:x11(2)),1:x12(2)-2,x12(3:x12(2))),title('牛顿迭代法');

x13=Secant(-1,1,1e-10);

x14=Secant(-4,5,1e-10);

x15=Secant(0,7,1e-10);

x16=Secant(-8,2,1e-10);

x17=StephensonIterative(-1,1e-10);

x18=StephensonIterative(-4,1e-10);

x19=StephensonIterative(4,1e-10);

x20=StephensonIterative(6,1e-10);

figure,

subplot(1,2,1),

plot(1:x3(2)-2,x3(3:x3(2)),1:x13(2)-2,x13(3:x13(2)),1:x14(2)-2,x14(3: x14(2)),1:x15(2)-2,x15(3:x15(2)),1:x16(2)-2,x16(3:x16(2))),title('牛顿割线法');

subplot(1,2,2),

plot(1:x4(2)-2,x4(3:x4(2)),1:x17(2)-2,x17(3:x17(2)),1:x18(2)-2,x18(3: x18(2)),1:x19(2)-2,x19(3:x19(2)),1:x20(2)-2,x20(3:x20(2))),title('史

蒂芬森迭代法');

实验结果：

1. 各个迭代值分布

图 1.1 不同迭代法下的得到的迭代值

迭代值的情况如下：

二分法牛顿迭代法

牛顿割线法

史蒂芬森迭代法

0 0 0 0 1.5000000000 0.2000000000 2.0000000000 1.3555555556 0.7500000000 0.3704918032 0.3333333333 0.9816165283 1.1250000000 0.5076442076 0.3807196801 0.9999460003 0.9375000000 0.6146189447 0.4982833419 0.9999999995

1.0312500000 0.6973869098 0.5704996333 0.9843750000 0.7615538091 0.6393806244 1.0078125000 0.8115411186 0.6942785879 0.9960937500 0.8506763857 0.7411692653 1.0019531250 0.8814482123 0.7802715997 0.9990234375

0.9057297400

0.8132927871

当二分法的初始区间选为[0,3]，误差限为1×10?10，牛顿迭代法初值选为0，

010203040

0.5

1.5二分法

050

100

0.5

牛顿迭代法

100

150

00.511.5

2牛顿割线法

02468

0.51

1.5史蒂芬森迭代法

误差限为1×10?10，牛顿割线法初始点为{0,2}，误差限为1×10?5，史蒂芬森迭代法初始点选为0，误差限为1×10?10，迭代情况如图所示。迭代次数分别为38次，100次，140次，9次。故而，史蒂芬森迭代法速度最快，效果最好。 2. 收敛情况

图 1.2 不同迭代法下迭代值得收敛情况

二分法收敛效果较差，牛顿迭代法和牛顿割线法相近，史蒂芬森迭代法收敛次数高于1，效果最好 3. 不同初值的收敛情况

-10

-5

二分法

100

0.40.5

0.6

0.70.8牛顿迭代法

-2

-1

1牛顿割线法

-0.4

-0.2

0.20.4史蒂芬森迭代法

图 1.3 二分法，牛顿迭代法下不同初值的收敛情况

图 1.4 牛顿割线法，史蒂芬森迭代法下不同初值的收敛情况

-4-3

-2

-1

二分法

50100150

-4-3

-2-10

345

6牛顿迭代法

020406080-8

-6-4

-20246

8牛顿割线法

02468

00.5

11.5

2.5

史蒂芬森迭代法

1.二分法的五个初始区间分别为[0,3],[?1,1],[?2,3],[0,4],[?4,4]；

2.牛顿迭代法的五个初始值分别为0,?2,?4,4,6；

3.牛顿割线法的五个初始区间分别为[0,2],[?1,1],[?4,5],[0,7],[?8,2]；

4.史蒂芬森迭代法的五个初始值分别为0,?1,?4,4,6；

由图可知，它们最终均达到收敛。

收敛性分析及结论：

1.二分法收敛较慢且不能求解崇根，但算法简单；此处牛顿法具有了平方收敛；

从迭代次数上看，牛顿割线法较牛顿法的多，所以收敛性较差，是超线性收敛；史蒂芬森迭代法收敛效果最好。

2.因为牛顿迭代法是局部的二次收敛，所以要注重初值的选取，本次实验中选

择的初值均得到了收敛，效果比较好。牛顿割线法也应注意初值的选取。

（第三章）

实验题目：

1.区间[?1,1]作等距划分：

x k=?1+k? ,?=2

,k=0,1,…,n

以x k为结点对函数f(x)=1

5+x2

进行插值逼近。

（1）分别取?=1,5,10,20,25用牛顿插值对f(x)进行逼近，并在同一坐标系下做出函数的图形，进行比较。写出插值函数对f(x)的逼近程度与节点个数的

关系，并分析原因；

（2）试用三次样条插值对f(x)进行逼近，在同一坐标下画出图形，观察样条插值函数对f(x)的逼近程度与节点个数的关系；

（3）整体插值有何局限性？如何避免？

2.已知一组数据如下，求其拟合曲线.

表 2.1数据表

i0 1 2 3 4 5 6 7 8 9 10 x i 2 3 4 7 8 10 11 14 16 18 19 y i106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2 （1）求以上数据形如 y(x)=c0+c1x+c2x2的拟合曲线及其平方误差；

（2）求以上数据形如 y(x)=ae?b

x的拟合曲线及其平方误差；

（3）通过观察（1）（2）的结果，写出你对数据拟合的认识. 问题分析：

题目除上述要求之外还有以下几点：

1.明确整体插值和分段插值的不同。牛顿插值多项式属于整体插值，三次样条

插值属于低次分段插值。

2.将结果在同一坐标下绘制出。但是为了方便分析节点个数对于插值效果的影

响，也可以单独绘制。

3.第二题中为了确定各个参数的大小，可以进行适当变换，转化为线性，运用

最小二乘法，得到拟合。

实验原理：

牛顿插值多项式：对于给定的插值节点a≤x0

N n(x)=a0+a1(x?x0)+a2(x?x0)(x?x1)+?

+a n(x?x0)(x?x1)…(x?x n?1)

使其满足插值条件

P n(x i)=y i=f(x i) i=0,1,…,n

这样得到的插值多项式 N n(x)称为 Newton 插值多项式。系数 a i为差商，可以通过构造差商表得到。

三次样条插值：三次样条插值函数 S(x)在每个小区间[x k?1,x k]上为三次多项式；在全进(a,b)上存在二阶连续导数；其次符合插值条件。Matlab 中存在内置的三次样条插值函数，命令为 spline .

实验内容：

第一题：

1.牛顿插值函数的构造

代码如下：

function f=Newton(x0,y0)

%牛顿多项式插值函数

syms x;

SZ=size(x0,2);

a(1)=y0(1);

y(:,1)=y0';

for j=2:SZ

nx1=1;

for i=1:SZ-j+1;

nx2=nx1+j-1;

y(i,j)=(y(i,j-1)-y(i+1,j-1))/(x0(nx1)-x0(nx2));

nx1=nx1+1;

end

f=y(1,1);

for j=2:SZ

ff=y(1,j);

for i=1:j-1

ff=ff*(x-x0(i));

end

f=f+ff;

end

2.牛顿和三次样条插值情况及比较：

代码如下：

clear all;

clc;

syms x;

fx=1/(5+x^2);

%牛顿多项式插值

x0=-1:0.01:1;

y0=subs(fx,x0);

n1=[1,5,10,20,25];

n2=[5,55,60,62,67];

h1=2./n1;

h2=2./n2;

x1=-1:h1(1):1; y1=subs(fx,x1); f1=Newton(x1,y1); y01=subs(f1,x0); x2=-1:h1(2):1; y2=subs(fx,x2); f2=Newton(x2,y2); y02=subs(f2,x0); x3=-1:h1(3):1; y3=subs(fx,x3); f3=Newton(x3,y3); y03=subs(f3,x0); x4=-1:h1(4):1; y4=subs(fx,x4); f4=Newton(x4,y4); y04=subs(f4,x0); x5=-1:h1(5):1; y5=subs(fx,x5); f5=Newton(x5,y5); y05=subs(f5,x0); figure,

plot(x0,y0,x0,y01,x0,y02,x0,y03,x0,y04,x0,y05),title('所有结果');

x6=-1:h2(1):1; y6=subs(fx,x6); f6=Newton(x6,y6); y06=subs(f6,x0); x7=-1:h2(2):1; y7=subs(fx,x7); f7=Newton(x7,y7); y07=subs(f7,x0); x8=-1:h2(3):1; y8=subs(fx,x8); f8=Newton(x8,y8); y08=subs(f8,x0);

x9=-1:h2(4):1; y9=subs(fx,x9); f9=Newton(x9,y9); y09=subs(f9,x0);

x10=-1:h2(5):1; y10=subs(fx,x10); f10=Newton(x10,y10);

y010=subs(f10,x0);

figure,

plot(x0,y0,x0,y06,x0,y07,x0,y08,x0,y09,x0,y010),title('龙格现象');

%三次样条插值spline自带命令

x0=-5:0.01:5;

y0=subs(fx,x0);

figure,

subplot(2,3,1),plot(x0,y0),title('精确图');

y11=spline(x1,y1,x0);subplot(2,3,2),plot(x0,y11),title('n=1');

y12=spline(x2,y2,x0);subplot(2,3,3),plot(x0,y12),title('n=5');

y13=spline(x3,y3,x0);subplot(2,3,4),plot(x0,y13),title('n=10');

y14=spline(x4,y4,x0);subplot(2,3,5),plot(x0,y14),title('n=20');

y15=spline(x5,y5,x0);subplot(2,3,6),plot(x0,y15),title('n=25'); figure,

plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15),title('所有结果');

第二题：

代码如下：

clear all;clc;

x=[2 3 4 7 8 10 11 14 16 18 19];

y=[106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2]; %Agenda 1

A=[(x.*x)' x' ones(11,1)];

A1=A'*A;

y1=A'*y';

c1=A1\y1

x1=2:0.1:19;

y1=polyval(c1,x1);

plot(x,y,'k*',x1,y1,'r'),gtext('曲线一'),hold on

y01=polyval(c1,x);

delta1_2=sum((y-y01).^2)

%Agenda 2

xx=-1./x;

yy=log(y);

B=[ones(11,1) xx'];

B1=B'*B;

yy1=B'*yy';

c2=B1\yy1;

syms s;

a=exp(c2(1));b=c2(2);

[a b]

f=a*exp(-b/s);

x2=x1;

y2=subs(f,x2);

plot(x2,y2,'b'),gtext('曲线二'); y02=subs(f,x);

delta2_2=sum((y-y02).^2)

实验结果：

第一题：

1.牛顿插值结果

图 2.1 不同节点数的牛顿插值多项式综合图

图 2.2 龙格现象

-1

-0.8-0.6-0.4-0.200.20.40.60.81

所有结果

-1-0.8-0.6-0.4-0.200.20.40.60.81

0.15

0.16

0.17

0.18

0.19

0.2

龙格现象

由图可知，

2. 三次样条插值结果

图 2.3 不同节点下的三次样条插值结果

-5

00.050.10.15

0.2精确图

-5

-101

2n=1

-5

00.2

0.40.6

0.8n=5

-5

500.20.40.6

0.8n=10

-5

00.20.40.6

0.8n=20

-5

00.20.40.6

0.8n=25

图 2.4 不同节点数的三次样条的综合图

由图可知，牛顿插值多项式在n 较小的时候如{5,10,20,25}差值效果良好，当n 变大时如{60,62,65,67,100,200}就出现了龙格现象，三次样条在各个子区间内为三次多项式，拟合效果好. 第二题

1. 第一问的多项式拟合得到的拟合曲线为

y (x )=106.2927+0.6264x ?0.0205x 2

平方误差为

δ1=2.7796

第二问的拟合曲线为

y (x )=111.4940e

0.0903x

平方误差为

δ2=0.4719

2. 拟合曲线如图所示

-5

-4-3-2-1012345

00.1

0.2

0.3

0.4

0.5

0.6

0.7

所有结果

2468101214161820

图 2.5拟合情况

从图中可以看出曲线二大体符合黑点的分布情况，拟合效果较好。

结论：

整体插值随着节点个数的增多，多项式的次数也在升高。高次多项式的插值会出现龙格现象，震荡剧烈，逼近效果不理想。但是当节点很多时，低次插值的精度又不够，所以为了避免这一局限性采用分段低次插值。其中三次样条插值有良好的收敛性和光滑性，效果较好。

数据拟合时，可以选择趋势相似的函数形式，在求解相关参量时，将函数进行适当变换，从而运用最小二乘法得到拟合结果。

（第四章）

实验题目：

设计区间分半求积算法、龙贝格求积算法和自适应辛普森求积算法的程序，观察n =1,10,100,500时，积分∫x 2cos (nx)dx 1

?1的结果，并给出相应的评价.

问题分析：

1. 分半求积算法即为复合梯度求积。

2. 因为n 越大，f (x )=x 2cos (nx)在[?1,1]内震荡地越严重，自适应辛普森求积算法很难适用，计算复杂度会很高。所以选取辛普森求积算法代替。

实验原理：

分半求积算法：将区间[a,b ]等分为n 个小区间，分点为

x k =a +k? ?=b ?a

,k =0,1,…,n

利用定积分性质，有

∫f (x )dx =∑∫

f(x)dx x k+1

x k

n?1

k=0b

每个小区间上均用梯形求积公式，有

∫f (x )dx =?2∑[f(x k )+f(x k+1)]+∑[??2

12f(ηk )]n?1

k=0

n?1

k=0

令

T n =?2∑[f(x k )+f(x k+1)]n?1

k=0

[f(a)+f(b)]+?∑f(x k )n?1

k=0

即为复合梯形公式。

龙贝格求积算法：将步长为?的复合梯形公式进行査尔逊外推，就得到龙贝格求积公式。

自适应辛普森求积算法：按照被奇函数在区间上的变化来安排求积结点的辛普森算法。

实验内容：

1.写得各个求积算法的函数。

代码如下：

（1）分半求积算法：

function

f=CompositeTrapezoida(n)

syms x;

y=x*x*cos(n*x);

m=1:100;

a=-1;

b=1;

for i=1:100

h=(b-a)/m(i);

k=0:m(i);

xk=a+k*h;

yk=subs(y,xk);

f(i)=0;

for j=1:m(i)

f(i)=f(i)+(yk(j)+yk(j+1));

end

f(i)=h/2*f(i);

end

（2）辛普森算法

function f=CompositeSimpson(n) syms x;

y=x*x*cos(n*x);

m=1:100;

a=-1;

b=1;

for i=1:100

h=(b-a)/(2*m(i));

k=0:(2*m(i)); xk=a+k*h;

yk=subs(y,xk);

f(i)=0;

for

j=1:m(i)f(i)=f(i)+(yk(2*j-1)+4*y k(2*j)+yk(2*j+1));

end

f(i)=h/3*f(i);

end

[m;f]';

end

（4）龙贝格求积算法

function f=Romberg(n)

epsilon=0.0001;

syms x;

y=x*x*cos(n*x);

a=-1;

b=1;

h=b-a;

fa=subs(y,a);

fb=subs(y,b);

T(1)=h/2*(fa+fb);

m=1;

while 1

h=h/2;

k=1:(2^m-1);

xk=a+k*h;

yk=subs(y,xk);

su=sum(yk);

S(1)=h/2*(fa+fb)+h*su;

for

i=1:mS(i+1)=S(i)+(S(i)-T(i))/(4^ i-1);

end

if(abs(S(m+1)-T(m))

end

T=S;

m=m+1;

end

f=S(m+1);

end

（5）自适应辛普森求积算法

function f=AdaptiveSimpson(n) epsilon=0.001;

syms x;

y=x*x*cos(n*x);

a=-1;

b=1;

h=b-a;

p=[a b];

p0=p;

ep=[epsilon];

m=0;

q=0;

f=0;

while 1

n1=length(ep);

n=length(p0);

if(n==1) break; end h=p0(2)-p0(1);

k=0:4;

yk=subs(y,p0(1)+k*h/4);

s0=h/6*(yk(1)+4*yk(3)+yk(5)); s1=h/12*(yk(1)+4*yk(2)+yk(3));

s2=h/12*(yk(3)+4*yk(4)+yk(5));

if(abs(s0-s1-s2)<15*ep(1))

f=f+s1+s2;

p0=p0(2:n);

if n1>=2

ep=ep(2:n1);

end

q=q+1;

else

m=m+1;

p0=[yk(1),yk(3),p0(2:n)]; if n1==1

ep=[ep(1)/2 ep(1)/2]; else

ep=[ep(1)/2 ep(1)/2

ep(2:n1)];

end

if q==0

p=p0;

else

p=[p(1:q) p0];

end

2.各个求积公式下的积分结果

代码如下:

clear all;clc;%复合梯形法

y11=CompositeTrapezoida(1);

y12=CompositeTrapezoida(10); y13=CompositeTrapezoida(100); y14=CompositeTrapezoida(500); [1:100;y11;y12;y13;y14]'

%复合辛普森法

y21=CompositeSimpson(1);

y22=CompositeSimpson(10);

y23=CompositeSimpson(100);

y24=CompositeSimpson(500); [2:2:200;y21;y22;y23;y24]'

%龙贝格法y31=Romberg(1);

y32=Romberg(10);

y33=Romberg(100);

y34=Romberg(500);

[y31 y32 y33 y34]

%自适应辛普森法

y41=AdaptiveSimpson(1)

%以下复杂度太高不出结果

%y42=AdaptiveSimpson(10); %y43=AdaptiveSimpson(100); %y44=AdaptiveSimpson(500); %[y41 y42 y43 y44]

实验结果：

表 3.1不同积分法得到的结果

结果评价：

龙贝格法的精度高，最接近精确解。因为被积分函数为f(x)=x2cos (nx)随着n的增大，震荡越剧烈。所以积分面积越来越小。

（第五章）

实验题目：

1.分别取步长 ?=0.5,0.1,0.05,0.01,… 用龙格-库塔法求解y′=ty 1

3 ,y(1)=1 ,

计算到t=2 , 并与精确解y(t)=[(t 2+2)

]

2比较，观察龙格-库塔法的收敛性.

2.分别取步长?=0.1,0.05,0.01,0.005,…用显示欧拉法和隐式欧拉法求解y′=

?50y ,y(0)=100 , 由结果分析算法的稳定性.

3.选择某常微分方程初值问题的数值方法计算∫sin t

t dt

的近似值，并保证有4位有效数字.