计算方法实验报告

计算方法实验报告（四）

（一）线性方程的迭代解法

一、实验问题

利用简单迭代法，两种加速技术，牛顿法，改进牛顿法，弦割法求解习题5-1，5-2，5-3中的一题，并尽可能准确。

选取5-3：求x 3−x 2−1=0在x=1.5附近的根。

二、问题的分析（描述算法的步骤等）

（1）简单迭代法算法：

给定初始近似值p 0,求p =φ(p )的解。

Step 1 令i=0；

Step 2 令p i+1=φ(p i )(计算p i+1)；

Step 3 如果p i+1=p i ,则迭代终止，否则重复Step 2。

（2）Aitken 加速法算法

Step 1 令k=0，利用简单迭代算法x k+1=φ(x k )得到迭代序列x k ；

Step 2 令x k ∗=x k -(x k −x k−1)2

x k −2x k−1+x k−2(计算x k ∗得到一个新的序列x k ∗，其中k=0,1,2…)；

Step 3 如果x k+1∗=x k ∗,则迭代终止，否则重复Step 2。

(3)插值加速法算法

Step 1 令k=0，利用简单迭代算法x k+1=φ(x k )得到迭代序列x k ；

Step 2 令x k ∗=x k +(x k −x k−1)(x k −x k+1)

x k−1−2x k +x k+1(计算x k ∗得到一个新的序列x k ∗，其中k=1,2

,3…)；

Step 3 如果x k+1∗=x k ∗,则迭代终止，否则重复Step 2。

(4)牛顿法算法

Step 1给定初始近似值x 0；

Step 2令x k+1=x k −f (x k

)f (x k ),其中k ∈N,计算得到x k 的 序列；

Step 3如果x k+1=x k ,则迭代终止，否则重复Step 2。

（5）改进牛顿法的算法

Step 1给定初始近似值x 0；

Step 2令x k+1=x k −k ′(k )′(k )√′2(k )"(k )(k ) ,其中

k ∈N,

迭代计算得到x k 的 序列；

Step 3如果x k+1=x k ,则迭代终止，否则重复Step 2。

（6）弦割法算法（双点弦割法）

Step 1给定初始近似值x 0，x 1；

Step 2令x k+1=x k −f (x k )(x k −x k−1)

f (x k )−f (x k−1)其中k ∈N,计算得到x k 的 序列；

Step 3如果x k+1=x k ,则迭代终止，否则重复Step 2。

三、程序设计

(1)简单迭代法

利用迭代公式x =√1+x 23进行迭代运算。

#include

#include

#include

double fun(double x)

{

double c=1+x*x;

return pow(c,1/3.0);

}

void main()

{

double x=1.5;

double y=0;

double D=1;

double e=0.001;

while(D>e)

{

D=0;

y=fun(x);

if(fabs(y-x)>=D)

{

D=fabs(y-x);

}

x=y;

}

cout<

}

(2) ）Aitken加速法源程序如下：

x1=1.5;

eps=0.0001;

y1=(1+x1^2)^(1/3);

z1=(1+y1^2)^(1/3);

x=z1-(z1-y1)^2/(z1-2*y1+x1);

while eps

x=x2;

x1=y1;

y1=(1+x1^2)^(1/3);

z1=(1+y1^2)^(1/3);

x2=z1-(z1-y1)^2/(z1-2*y1+x1); n=n+1;

end

fprintf('迭代次数 n=%.0f\n',n);

fprintf('x2=%.5f\n',x2)

（3）插值加速法源程序如下：

x1=0;

x2=1.5;

eps=0.0000001;

y1=0;

z1=0;

n=0;

while eps

x1=x2;

y1=(1+x1^2)^(1/3);

z1=(1+y1^2)^(1/3);

x2=z1+(z1-y1)^2/(z1-2*y1+x1); n=n+1;

end

fprintf('迭代次数 n=%.0f\n',n);

fprintf('x2=%.5f\n',x2)

（4）牛顿法：

进行迭代运算

利用公式x k+1=x k−f(x k)

f(x k)

程序设计如下：

#include

#include

double fun(double x)

{

double a=2*pow(x,3.0)-pow(x,2.0)+1;

double b=3*pow(x,2.0)-2*x;

return a/b;

}

void main()

{

double x=1.5;

double y=0;

double D=1;

double e=0.001;

double f=0;

while(D>e)

{

D=0;

y=fun(x);

if(fabs(y-x)>=D)

{

D=fabs(y-x);

}

x=y;

f++;

}

cout<