数值分析上机实践

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

1、分别用牛顿法，及牛顿-Steffensen法

(1) 求ln(x+sinx)=0的根。初值x0分别取0.1, 1,1.5, 2, 4进行计算。 (2) 求sinx=0的根。初值x0分别取1，1.4，1.6, 1.8，3进行计算。 分析其中遇到的现象与问题。

（1）、newton法

(1)、对应M文件为newton1如下： function x=newton1(x0)

x=x0-((x0+sin(x0))*log(x0+sin(x0)))/(1+cos(x0)); while (x-x0>0)|(x-x0<0) x0=x;

x=x0-((x0+sin(x0))*log(x0+sin(x0)))/(1+cos(x0)); end

>> newton1(0.1)

ans =

0.5110

>> newton1(1)

ans =

0.5110

1.5、2.4作为初值不适合，计算结果出现复数。

(2)对应M文件为newton2如下： function x=newton2(x0) x=x0-sin(x0)/cos(x0); while (x-x0>0)|(x-x0<0) x0=x;

x=x0-sin(x0)/cos(x0); end

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

>> newton2(1)

ans =

>> newton2(1.4)

ans =

3.1416

>> newton2(1.6)

ans =

31.4159

>> newton2(1.8)

ans =

6.2832

>> newton2(3)

ans =

3.1416

>> newton2(3.3)

ans =

3.1416

>> newton2(3.9)

ans =

3.1416

初始值取值不同，有可能造成迭代公式发散。计算不出正确值，或是计算得出错误的值。

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

（2）、newton-steffensen加速法

(1)、对应的M文件为NS111如下：

function x=NS111(x0)

y0=x0-((x0+sin(x0))*log(x0+sin(x0)))/(1+cos(x0)); z0=y0-((y0+sin(y0))*log(y0+sin(y0)))/(1+cos(y0)); x=x0-(y0-x0)*(y0-z0)/(z0-2*y0+x0); while (abs(x-x0)>0.00001) x0=x;

y0=x0-((x0+sin(x0))*log(x0+sin(x0)))/(1+cos(x0)); z0=y0-((y0+sin(y0))*log(y0+sin(y0)))/(1+cos(y0)); x=x0-(y0-x0)*(y0-z0)/(z0-2*y0+x0); end

x=x0

>> NS111(0.1) x =

-10.7295 - 3.1908i

ans =

-10.7295 - 3.1908i

>> NS111(1) x =

0.5140

ans =

0.5140

>> NS111(1.5) x =

0.5140 - 0.0000i

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

ans =

0.5140 - 0.0000i

>> NS111(2.4) x =

2.9105e+002 -4.0048e-001i

ans =

2.9105e+002 -4.0048e-001i

（2）、对应M文件为NS2如下： function x=NS2(x0) y0=x0-sin(x0)/cos(x0); z0=y0-sin(y0)/cos(y0);

x=x0-(y0-x0)*(y0-z0)/(z0-2*y0+x0); while (abs(x-x0)>0.00001) x0=x;

y0=x0-sin(x0)/cos(x0); z0=y0-sin(y0)/cos(y0);

x=x0-(y0-x0)*(y0-z0)/(z0-2*y0+x0); end

x=x0

>> NS2(1) x =

0.0311

ans =

0.0311

>> NS2(1.4)

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

x =

-0.0311

ans =

-0.0311

>> NS2(1.6) x =

3.1727

ans =

3.1727

>> NS2(1.8) x =

9.3937

ans =

9.3937

>> NS2(3) x =

3.1105

ans =

3.1105 >> NS2(3.3)

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

3.1727

ans =

3.1727

两种算法算出的X值有略微的差别，说明两者的计算精度不同.

2 松弛因子对SOR法收敛速度的影响。

用SOR法求解方程组Ax=b,其中

要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代，记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w<=0或w>=2会有什么影响？

解：本题的主要目的为分析松弛因子对收敛速度的影响,在编写程序是需要输出迭代次数,

程序能够接受不同的松弛因子w,而且可以算出方程的解.

本题方程组很特殊,无论是几阶,其解都是1. 而且无论是几阶,都可以将SOR迭代公

式写成下式 : (

)

使用上式不难编出程序,本程序可以输入方程阶数和w值,可以求出迭代次数和解向量,程序的初始值设置为0向量,程序中设置迭代终止条件||x(k+1)-x(k)||1<10-6 ,本程序只适于上式,不具普遍性.由于本题方程组很特殊,所以程序既没有存系数矩阵,也

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

没有用向量存系数.源代码见附录A .

下面分析5到9阶的方程,w取为1.1到1.9,使用程序观察迭代次数.结果总结在

表8中.

表8 迭代次数计算汇总

从上面的结果可以看出,对于本方程组,随着w的增大,大大地降低了迭代速度,方程阶数的提高也有增大迭代次数的趋势.

下面分析w<=0 的情况.还是取5到9阶的方程,w分别取为-0.1 -0.5 -1.0 -2.0 四个值,结果均不收敛.

下面分析w>=2 的情况. 还是取5到9阶的方程,w分别取为2.0 2.1 2.5 3 四个值,结果均不收敛.

从书上的理论可知,SOR迭代法收敛的一个必要条件是0<w<2 ,程序验算的结果证明了这个理论.

总结与分析

1. SOR迭代法的关键还是写出收敛的迭代格式,C++程序只负责迭代计算. 2. 从本题来看w取为1.1时迭代速度最快,取为1.9时收敛速度最慢.

3. SOR方法是求解大型稀疏矩阵方程组的有效方法之一,如果遇到系数矩阵是稀疏

矩阵,那么应该想到这种方法.

数值分析一门课是理论与实践相结合的一门课,只有理论和实践相结合才能达到好的效果,本文给出了几个上机的实例,希望能给初学者以好的帮助

4. 对本题这种三对角方程组,追赶法更加适宜.

SOR迭代法程序（C++） #include<iostream> #include<cmath> #include<iomanip> #include<vector> using namespace std;

double norm1(const vector<double>& v1,const vector<double>& v2) { int i; double ans=0.0; for(i=0;i<=v1.size()-1;i++) ans+=fabs(v1[i]-v2[i]); return ans; }

int main() { int i,n,js=0; double w; cin>>n; cin>>w; vector<double> x1(n),x2(n); fill(x1.begin(),x1.end(),0.0); fill(x2.begin(),x2.end(),0.0); while(1) { x2[0]+=(-(w/4)*(-3-x2[1]+4*x2[0])); for(i=1;i<=x2.size()-2;i++) x2[i]+=(-(w/4)*(-2-x2[i-1]+4*x2[i]-x2[i+1])); x2[x2.size()-1]+=(-(w/4)*(-3+4*x2[x2.size()-1]-x2[x2.size()-2])); js++; if(norm1(x1,x2)<1e-6) break; else x1=x2; } cout<<"迭代次数:"<<js<<endl; cout<<"解向量:"<<endl; for(i=0;i<=x2.size()-1;i++) cout<<fixed<<setprecision(8)<<x2[i]<<endl; return 0;