数值分析实验报告

数值分析实验报告

《数值分析》实验报告

班级：

学号： 姓名：

1

数值分析实验报告

课题1 解线性方程组的直接算法

一、问题提出

给出下列几个不同类型的线性方程组，请用适当算法计算其解。

1、设线性方程组

2?3?1210000??x1??4?5??8??x??12?6?5?36501002???????4?3?2?2?132?1031??x3???????x0?215?13?11942???4?????42?3?6?167?3323??x5?????=??

6?8571726?35??x6??8?46??0?13?2?13?425301??x7???????2?122??x8??1610?11?91734?38????4?19?62?713920124??x9???????0?18?3?24?863?1??0????21???x10??x*= -1, 0, 1, 2, 0, 3, 1, -1, 2 )T 2、设对称正定阵系数阵线方程组

2?402400??x1??0??4?2??x???6?2?1?21320???2?????4?1141?8?356??x3??20???????x0?216?1?4?3323???4?= ??

?21?8?1224?10?3??x5??9???????43?3?44111?4?22x???6????025?3?101142??x7???15???????063?3?4219???0???45???x8???x* = ( 1, -1, 0, 2, 1, -1, 0, 2 )T

3、三对角形线性方程组

?4?100000000??x1??7???14?10000000??x??5?2???????0?14?1000000??x3???13???????x00?14?1000002???4????000?14?10000??x5??6? ???? = ??

0000?14?1000?12x???6????00000?14?100??x??14?7???????000000?14?10??x8???4??0000000?14?1??x??5????9??????00000000?14?????5???x10??x*= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )T

2

数值分析实验报告

二、要求

1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）； 2、编出算法通用程序

3、在应用Gauss消去法时，尽可能利用相应程序输出系数矩阵的三角分解式。

三、目的和意义

1、通过该课题的实验，掌握模块化结构程序设计方法的优点； 2、掌握解各类线性方程组的直接算法，了解各种方法的特点； 3、体会高斯消去法选主元的必要性。； 四、流程图

3

入口 定义各种变量 请输入方程维数n 请输入方程系数矩阵A 请输入矩阵B 请输入精度要求 \请输入最大迭代次数 请输入初始向量 数值分析实验报告

用高斯消去法解方程组 用平方根法解方程用追赶法解方程结束

五、源程序代码

1、Gauss顺序消去法 #include \void main() {

double a[10][10]={4,2,-3,-1,2,1,0,0,0,0, 8,6,-5,-3,6,5,0,1,0,0, 4,2,-2,-1,3,2,-1,0,3,1, 0,-2,1,5,-1,3,-1,1,9,4, -4,2,6,-1,6,7,-3,3,2,3, 8,6,-8,5,7,17,2,6,-3,5, 0,2,-1,3,-4,2,5,3,0,1, 16,10,-11,-9,17,34,2,-1,2,2, 4,6,2,-7,13,9,2,0,12,4, 0,0,-1,8,-3,-24,-8,6,3,-1}; double b[10]={5,12,3,2,3,46,13,38,19,-21}; double s; int i,j,k;

for(i=0;i<10;i++) {

for (k=i+1;k<10;k++) {

s=-a[k][i]/a[i][i]; for(j=0;j<10;j++) {

a[k][j]=s*a[i][j]+a[k][j]; }

b[k]=b[i]*s+b[k];

4

数值分析实验报告

} }

for(i=9;i>=0;i--) { s=0;

for(j=9;j>i;j--) s=s+a[i][j]*b[j]; b[i]=(b[i]-s)/a[i][i]; }

for(i=0;i<=9;i++)

printf(\}

2、平方根法 #include \ #include \ void main() {

FILE *f;

double a[8][8],l[8][8]; double s; int i,j,k,n;

f=fopen(\ fscanf(f,\ for(i=1;i<=n;i++) for(j=1;j<=i;j++)

fscanf(f,\ fclose(f);

for(i=1;i<=n;i++) {

for(j=1;j<=i-1;j++) {

for(s=0,k=1;k<=j-1;k++) s+=l[i][k]*l[k][j]; l[i][j]=(a[i][j]-s)/l[j][j]; } s=0;

for(s=0,k=1; k<=i-1; k++) s+=l[i][k]*l[i][k]; l[i][i]=sqrt(a[i][i]-s); }

for(i=1;i<=n;i++) {

5

数值分析实验报告

for(j=1;j<=i;j++) printf(\ l[i][j]); } }

sqrt.dat数据文件格式如下： 8 4 2 2

-4 -1 14

0 -2 1 6

2 1 -8 -1 22

4 3 -3 -4 4 11

0 2 5 -3 -10 1 14 0 0 6 3 -3 -4 2 19

3、追赶法

#include \ void main() {

FILE *f;

double a[10],b[10],c[10],d[10]; double t; int i,n;

f=fopen(\ fscanf(f,\

fscanf(f,\ for(i=1;i<=n-1;i++) {

fscanf(f,\ }

fscanf(f,\ fclose(f);

c[0]=c[0]/b[0]; d[0]=d[0]/b[0];

for(i=1;i<=n-1;i++) {

t=b[i]-c[i-1]*a[i]; c[i]=c[i]/t;

d[i]=(d[i]-d[i-1]*a[i])/t; }

d[n]=(d[n]-d[n-1]*a[n])/(b[n]-c[n-1]*a[n]); for(i=n-1;i>=0;i--) d[i]=d[i]-c[i]*d[i+1];

6

数值分析实验报告

for(i=0;i<=n;i++)

printf(\ }

zgf.dat数据文件格式如下： 10

4 -1 7 -1 4 -1 5 -1 4 -1 -13 -1 4 -1 2 -1 4 -1 6 -1 4 -1 -12 -1 4 -1 14 -1 4 -1 -4 -1 4 -1 5 -1 4 -5 六、运行结果 第一个方程组：

第二个方程组

7

数值分析实验报告

第三个方程组：

课题3 线性方程组的迭代法

一、问题提出

对课题一所列的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidol迭代法和SOR方法计算其解。

8

数值分析实验报告

二、要求

1、应用迭代法求解线性方程组，并能与直接法做以比较；

2、分别对不同精度要求，如??10?3,10?4,10?5由迭代次数体会该迭代法的收敛快慢；

3、对方程组2，3使用SOR方法时，选取松弛因子?=0.8，0.9，1，1.1，1.2等，试观察对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者； 三、目的和意义

1、通过上机计算体会迭代法求解线性方程组的特点，掌握求解各类线性方程组的迭代法；

3、体会上机计算时，终止准则x(k?1)?x(k)?< ? 对控制迭代精度的有效性；

4、体会初始解 x(0)，松弛因子的选取，对迭代收敛速度的影响。 四、流程图

五、源程序代码

#include using namespace std;

#define N 40 const int n=10;

int jacobi(float *p,float b[],float X[],float x[],int n); int GS(float *p,float b[],float X[],float x[],int n); int SOR(float *p,float b[],float X[],float x[],int n);

9

数值分析实验报告

void print(float *a,int r);

void main() {

float A[10][10]={4,2,-3,-1,2,1,0,0,0,0, 8,6,-5,-3,6,5,0,1,0,0, 4,2,-2,-1,3,2,-1,0,3,1, 0,-2,1,5,-1,3,-1,1,9,4, -4,2,6,-1,6,7,-3,3,2,3, 8,6,-8,5,7,17,2,6,-3,5, 0,2,-1,3,-4,2,5,3,0,1, 16,10,-11,-9,17,34,2,-1,2,2, 4,6,2,-7,13,9,2,0,12,4, 0,0,-1,8,-3,-24,-8,6,3,-1}; float a[10]={5,12,3,2,3,46,13,38,19,-21}; float X1[10]={ 0,0,0,0,0,0,0,0,0,0}; float x1[10];

float B[8][8]={4,2,-4,0,2,4,0,0,

2,2,-1,-2,1,3,2,0, -4,-1,14,1,-8,-3,5,6,

0,-2,1,6,-1,-4,-3,3, 2,1,-8,-1,22,4,-10,-3, 4,3,-3,-4,4,11,1,-4, 0,2,5,-3,-10,1,14,2, 0,0,6,3,-3,-4,2,19}; float b[8]={0,-6,20,23,9,-22,-15,45}; float X2[8]={ 0,0,0,0,0,0,0,0}; float x2[8];

float C[10][10]={4,-1,0,0,0,0,0,0,0,0, -1,4,-1,0,0,0,0,0,0,0, 0,-1,4,-1,0,0,0,0,0,0, 0,0,-1,4,-1,0,0,0,0,0, 0,0,0,-1,4,-1,0,0,0,0, 0,0,0,0,-1,4,-1,0,0,0, 0,0,0,0,0,-1,4,-1,0,0, 0,0,0,0,0,0,-1,4,-1,0, 0,0,0,0,0,0,0,-1,4,-1, 0,0,0,0,0,0,0,0,-1,4}; float c[10]={7,5,-13,2,6,-12,14,-4,5,-5}; float x3[10];

float X3[10]={ 0,0,0,0,0,0,0,0,0,0}; float *p[3]; p[0]=&A[0][0]; p[1]=&B[0][0]; p[2]=&C[0][0];

cout<<\迭代法解第一个方程:\ jacobi(p[0],a,X1,x1,10);

10

数值分析实验报告

cout<<\迭代法解第二个方程:\ jacobi(p[1],b,X2,x2,8);

cout<<\迭代法解第三个方程:\ jacobi(p[2],c,X3,x3,10);

cout<<\迭代法解第一个方程:\ GS(p[0],a,X1,x1,10);

cout<<\迭代法解第二个方程:\ GS(p[1],b,X2,x2,8);

cout<<\迭代法解第三个方程:\ GS(p[2],c,X3,x3,10);

cout<<\迭代法解第一个方程:\ SOR(p[0],a,X1,x1,10);

cout<<\迭代法解第二个方程:\ SOR(p[1],b,X2,x2,8);

cout<<\迭代法解第三个方程:\ SOR(p[2],c,X3,x3,10); }

int jacobi(float *p,float b[],float X[],float x[],int n) {

int k,i,j; float m,R,r,e;

cout<<\请输入精度e:\ cin>>e;

for(k=0;k

R=0;

for(i=0;i

m=0;

for(j=0;j

m=m+(*(p+i*n+j))*X[j]; x[i]=X[i]+(b[i]-m)/(*(p+i*n+i)); r=x[i]-X[i]; if(r<0)

r=X[i]-x[i]; if(r>R) R=r; }

if(R<=e&&R>0) {

print(x,n);

cout<<\迭代次数为:\ return k; }

for(j=0;j<10;j++) X[j]=x[j]; }

11

// jacobi 数值分析实验报告

print(x,n);

cout<<\迭代次数为:\

cout<<\方程解发散,无法用Jacobi方法解此方程!\ return 0; }

int GS(float *p,float b[],float X[],float x[],int n) // GS { int i,j,k; float t,R,r,e;

cout<<\请输入精度e:\ cin>>e;

for(k=0;k

for(i=0;i

for(j=0;j

if(ji) t+=(*(p+i*n+j))*X[j]; }

x[i]=(b[i]-t)/(*(p+i*n+i)); }

for(i=0;i<10;i++) {

r=x[i]-X[i]; if(r<0)

r=X[i]-x[i]; if(r>R) R=r; }

if(R<=e&&R>0) {

print(x,n);

cout<<\迭代次数为:\ return k; }

for(j=0;j<8;j++) X[j]=x[j]; } print(x,n);

cout<<\迭代次数为:\

cout<<\方程解发散,无法用Gauss-Seidel方法解此方程!\ return 0; }

12

数值分析实验报告

int SOR(float *p,float b[],float X[],float x[],int n) // SOR {

int i,j,k;

float t,R,r,e,w;

cout<<\请输入松弛因子w(0>w;

cout<<\请输入精度e:\ cin>>e;

for(i=0;i

R=0;

for(i=0;i

t=0;

for(j=0;jR) R=r; }

if(R<=e&&R>0) {

print(x,n);

cout<<\迭代次数为:\ return k; } } print(x,n);

cout<<\迭代次数为:\

cout<<\方程解发散,无法用XOR方法解此方程!\ return 0; }

void print(float *a,int n) {

int j;

float *t=a; cout<<\ for(j=0;j

13

数值分析实验报告

六、运行结果

Jacobi迭代法求解截图

Gauss-Seidol迭代求解截图

SOR迭代法

14

数值分析实验报告

七、结果分析

由结果看出，解线性方程组的迭代方法与直接方法不同，即使在计算过程中无舍入误差，迭代法也难以获得精确解，迭代法是一类逐次近似的方法；分别取不同精度要求，??10?3,10?4,10?5，由迭代次数发现精度越高，迭代次数越多，收敛越慢；对方程组3使用SOR方法时，选取松弛因子?=0.8，0.9，1，1.1，1.2不同，对收敛快慢有影响，当松弛因子?=0.8时，收敛最快。

15

数值分析实验报告

课题6 函数插值方法

一、问题提出：

对于给定的函数y?f(x)的n+1个节点值yj?f(xj),j?0,1,?,n。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。

（1）给定数据如下：

0.4 0.55 xj 0.65 0.80 0.95 1.05 1.25382 yj

0.41075 0.57815 0.69675 0.90 1.00 求5次Lagrange多项式和分段2次插值多项式，计算f(0.596),f(0.99) 的值。

（2） 给定数据如下： 1 2 xj 3 4 5 6 7 yj 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，并计算的f(1.8)的值。

二、要求 ：

1、利用Lagrange插值公式编写出插值多项式程序；

2、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； 3、绘制插值多项式的曲线，观察其光滑性。 三、目的和意义：

1、 会用基本的插值方法求函数的近似表达式；

2、 了解插值多项式和分段插值多项式各自的优缺点； 3、 掌握插值方法的程序编制； 四、源程序代码： #include void main() {

printf(\计算机0906班******曲丽******20093062***********\\n\\n\

float x[]={0.4,0.55,0.65,0.80,0.95,1.05};

float y[]={0.41075,0.57815,0.69675,0.90,1.00,1.25382}; float l[7]; float p,h; int i,j,k;

float x1=0.596,x2=0.99;

16

数值分析实验报告

17

float y1,y2; float s[4];

for(i=0;i<6;i++) {

p=1;

for(j=0;j<6;j++)

if(j!=i) p=p*(x[i]-x[j]); l[i]=y[i]/p; }

l[6]=0;

printf(\对给定的6个节点构造5次Lagrange插值多项式为：\\n\printf(\for(i=0;i<6;i++) {

if(l[i]!=0) {

printf(\ for(j=0;j<6;j++)

if(j!=i) printf(\ if(l[i+1]>0) printf(\ } }

printf(\y1=0; y2=0;

for(i=0;i<6;i++) {

p=1;

for(j=0;j<6;j++)

if(j!=i) p=p*(x1-x[j]); y1+=l[i]*p; }

for(i=0;i<6;i++) {

p=1;

for(j=0;j<6;j++)

if(j!=i) p=p*(x2-x[j]); y2+=l[i]*p; }

printf(\printf(\h=x[5]-x[1];

s[0]=2*y[1]/(h*h); s[1]=-4*y[3]/(h*h);

数值分析实验报告

s[2]=2*y[5]/(h*h); s[3]=0;

printf(\对给定数据选取区间[0.55,1.05],构造分段二次Lagrange插值多项式为：\\n\ printf(\ for(i=0;i<3;i++) {

if(s[i]!=0) {

printf(\ for(j=0;j<3;j++)

if(j!=i) printf(\ if(s[i+1]>0) printf(\ } }

printf(\ y1=0; y2=0;

for(i=0;i<3;i++) {

p=1;

for(j=0;j<3;j++)

if(j!=i) p=p*(x1-x[2*j+1]); y1+=s[i]*p; }

for(i=0;i<3;i++) {

p=1;

for(j=0;j<3;j++)

if(j!=i) p=p*(x2-x[2*j+1]); y2+=s[i]*p; }

printf(\ printf(\ }

五、运行结果

18

六

数值分析实验报告

数据拟合曲线：

L(6)(x)

19

数值分析实验报告

L(5)(x)

S(2)(x)

20

数值分析实验报告

printf(\请输入积分下限b=\scanf(\

printf(\请输入允许的最大误差\ printf(\scanf(\

T[0][0]=(b-a)/2*(f(a)+f(b));

T[0][1]=0.5*T[0][0]+0.5*(b-a)*f((a+b)/2); T[1][0]=(4*T[0][1]-T[0][0])/(4-1); i=1;

while (fabs(T[i][0]-T[i-1][0])>e) { i++;

temp=pow(2,i-1); sum=0;

for(j=1;j<=temp;j++)

{ sum=sum+f(a+(2*j-1)*(b-a)/pow(2,i)); } T[0][i]=0.5*(T[0][i-1]+(b-a)/pow(2,i-1)*sum); for(m=1;m<=i;m++) { k=i-m;

T[m][k]=(pow(4,m)*T[m-1][k+1]-T[m-1][k])/(pow(4,m)-1); } } q=0; help=i; for(p=0;p<=i;p++) { printf(\for(q=0;q<=help;q++)

{ printf(\} help--; }

26

数值分析实验报告

T[0][0]=T[i][0];

printf(\为所求积分!\\n\}

五、运行结果 复合梯形求积

复合Simpson求积

27

数值分析实验报告

Romberg求积公式

六、结果分析

分别用梯形求积公式和复合Simpson求积公式对第二题求积分，比较结果可以看出，当等分的区间n为20时，复合Simpson求积公式的求解精度远远高于复合梯形求积公式，收敛速度十分快，Romberg求积公式为复合梯形积分的改进，收敛速度最快。通过不同的等分区间，可以看出取不同步长h?? b?a ?/n，求解精度不同，步长h越小，求解精度越高。给定精度要求?=0.000001，梯形求积公式n取到210时，才能达到所要求的精度，h取1/210,；复合Simpson求积公式n=10时，就达到了所要求的精度，h取0.1；Romberg求积公式n=8时，就达到了所要求的精度，h取1/8。

28

本文来源：https://www.bwwdw.com/article/xfh8.html