fortran数值计算基础

数值计算基础

目录

实验一直接法解线性方程组的 (2)

实验二插值方法 (11)

实验三数值积分 (5)

实验四常微分方程的数值解 (7)

实验五迭代法解线性方程组与非线性方程 (9)

实验一 直接法解线性方程组

一、实验目的

掌握全选主元消去法与高斯-塞德尔法解线性方程组。

二、实验内容

分别写出Guass 列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。

1、用Guass 列选主元消去法求解方程组

??????????=????????????????????--5.58.37.33.47.11.85.16.93.51.53.25.2321x x x

2、用追赶法求解方程组

??

?

??????

???????-=????????????????????????????????-----000010210000210000210000210000

254321x x x x x 三、实验仪器设备与材料

主流微型计算机

四、实验原理

1、Guass 列选主元消去法 对于AX =B

1）、消元过程：将（A|B ）进行变换为)~|~(B A ，其中A ~

是上三角矩阵。即：

????

??

?

??→

??????? ??n nn

n

n n nn

n n n

n b a b a b a a b a a a b a a a b a a a 0

01

01221112

2

1

222221111211 k 从1到n-1

a 、 列选主元

选取第k 列中绝对值最大元素ik n

i k a ≤≤max 作为主元。

b 、 换行

i

k ij kj b b n k j a a ?+=?,,1,

c 、 归一化 k

kk k kj kk kj b a b n k j a a a ?+=?/,,1,/

d 、 消元

n

k i b b a b n k j n k i a a a a i k ik i ij kj ik ij ,,1,,,1;,,1, +=?-+=+=?-

2）、回代过程：由)~

|~(B A 解出11,,,x x x n n -。

1

,2,,1,/1

-=?-

?∑+=n k x x a

b x a b k n

k j j kj

k n

nn n

2、追赶法 线性方程组为：

??

????????

? ??=??????????? ????????????? ?

?-----n n n n n n

n n n f f f f f x x x x x a b c a b c a b c a b c a 132********

3

3

3

22211

做LU 分解为：

??

??

???

?

?

?

??=??????????

?

?=-1111,1213

3221n n n R L βββαγαγαγα

分解公式：

????

?

????

-===-====-)1,,2,1(),,3,2(,),,3,2(1

11n i c n i b b n i a i i i i i i i i i αββγααγ 则

?

?

?==?=?=y Ux f

Ly f LUx f Ax 回代公式：

???

???

?

=-==-)

,,3,2(1

111n i y f y f y i i i i i αγα

??

?--=-==+)

1,,2,1(1

n n i x y x y x i i i i n n β

五、实验步骤

1、理解并掌握全选主元消去法与高斯-塞德尔迭代法公式；

2、画出全选主元消去法与高斯-塞德尔迭代法的流程图

3、使用C 语言编写出相应的程序并调试验证通过

六、实验报告要求

1、统一使用《武汉科技大学实验报告》本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。

2、源程序需打印后粘贴在实验报告册内；

3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。

七、实验注意事项

注意如何定义数据结构以保存矩阵和解以降低算法的复杂性。

八、思考题

若使用全主元消去法，在编程中应如何记录保存对于未知数的调换。

实验二 插值方法

一、实验目的

掌握拉格郎日插值法与牛顿插值法构造插值多项式。

二、实验内容

分别写出拉格郎日插值法与牛顿插值法的算法，编写程序上机调试出结果，要求所编程序适用于任何一组插值节点，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。

已知下列函数表

求x=0.5635时的函数值。三、实验仪器设备与材料

主流微型计算机

四、实验原理

已知n 个插值节点的函数值，则可由拉格郎日插值公式与牛顿插值公式构造出插值多项式，从而由该插值多项式求出所要求点的函数值。拉格郎日插值公式与牛顿插值公式如下：

1、Lagrange 插值公式

)

()(...)()()(0

1100x l y y x l y x l y x l x L n

k k k n n n ∑==+++=∏

≠=+-+---=----------=n k

j j j

k j n k k k k k k k n k k k x x x x x x x x x x x x x x x x x x x x x x x x x l 011101110)())(())(()

())(())(()( 2、Newton 插值公式

)

())(](,,[)

)(](,,[)](,[)()(11010102100100----++--+-+=n n n x x x x x x x x x f x x x x x x x f x x x x f x f x N

五、实验步骤

1、理解并掌握拉格郎日插值法与牛顿插值法的公式；

2、画出拉格郎日插值法与牛顿插值法算法的流程图；

3、使用C 语言编写出相应的程序并调试验证通过。

六、实验报告要求

1、统一使用《武汉科技大学实验报告》本书写，实验报告的内容要求有：实验目的、

实验内容、程序流程图、源程序、运行结果及实验小结六个部分。

2、源程序需打印后粘贴在实验报告册内；

3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。

七、实验注意事项

Newton插值法在编程时应注意定义何种数据结构以保存差商。

八、思考题

比较Lagrange插值法与Newton插值法的异同。

实验三 数值积分

一、实验目的

掌握复化梯形法与龙贝格法计算定积分。

二、实验内容

分别写出变步长梯形法与Romberge 法计算定积分的算法，编写程序上机调试出结果，要求所编程序适用于任何类型的定积分，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。

求

00001.0,sin 1

0≤?εdx x x

。

三、实验仪器设备与材料

主流微型计算机

四、实验原理

通过变步长梯形法与龙贝格法，我们只要知道已知n 个求积节点的函数值，则可由相应的公式求出该函数的积分值，从而不需要求该函数的原函数。变步长梯形法与龙贝格法公式如下：

1、变步长梯形法

∑∑-=-=+++=+=1

11

01)]()(2)([2)]

()([2

n i i n i i i n b f x f a f h

x f x f h

T

∑-=++=1

2/12)(221n i i n n

x f h T T

用ε≤-n n T T 2来控制精度 2、龙贝格法

∑-=++=1

2/12)(221n i i n n

x f h T T

n n n T T S 31342-=

n

n n S S C 15115162-= n

n n C C R 63163642-=

用ε≤-n n R R 2来控制精度

五、实验步骤

1、理解并掌握变步长梯形法与龙贝格法的公式；

2、画出变步长梯形法与龙贝格法的流程图

3、使用C 语言编写出相应的程序并调试验证通过

六、实验报告要求

1、统一使用《武汉科技大学实验报告》本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。

2、源程序需打印后粘贴在实验报告册内；

3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。

七、实验注意事项

在

?1

0sin dx x x

积分中，被积函数在x=0点函数值为1，对该点在程序设计中应注意对其

的定义。

八、思考题

使用复化梯形法与复化Simpson 法来计算该问题有何缺点？

实验四 常微分方程的数值解

一、实验目的

掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。

二、实验内容

分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。

求??

?≤≤=-=')

50(2

)0(2

x y xy y 步长h=0.25。

三、实验仪器设备与材料

主流微型计算机

四、实验原理

常微分方程的数值解主要采用“步进式”，即求解过程顺着节点排列次序一步一步向前推进，在单步法中改进欧拉法和四阶龙格-库塔法公式如下：

1、改进欧拉法

),(1n n n n y x hf y y +=+

)]

,(),([2111+++++=n n n n n n y x f y x f h

y y

2、四阶龙格-库塔法

????

?

??

???

???

++=++=++==++++=+)

,()2

,2()2,2(),()22(6

3423

12143211

hk y h x f k k h y h x f k k h y h x f k y x f k k k k k h y y n n n n n

n n n n n 五、实验步骤

1、理解并掌握改进欧拉法与四阶龙格-库塔法的公式；

2、画出改进欧拉法与四阶龙格-库塔法的流程图

3、使用C 语言编写出相应的程序并调试验证通过

六、实验报告要求

1、统一使用《武汉科技大学实验报告》本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。

2、源程序需打印后粘贴在实验报告册内；

3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。

七、实验注意事项

??

?≤≤=-=')

50(2

)0(2

x y xy y 的精确解为)1/(22x y +=，通过调整步长，观察结果的

精度的变化

八、思考题

如何对四阶龙格-库塔法进行改进，以保证结果的精度。

实验五 迭代法解线性方程组与非线性方程

一、实验目的

掌握高斯-塞德尔迭代法求解线性方程组与牛顿迭代法求方程根。

二、实验内容

分别写出高斯-塞德尔迭代法与牛顿迭代法的算法，编写程序上机调试出结果，要求所编程序适用于任何一个方程的求根，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。

1、高斯-塞德尔迭代法求解线性方程组

?????

???????-=????????????????????????----01741323151122231592127

4321x x x x 2、用牛顿迭代法求方程013=--x x 的近似根，00001.0≤ε，牛顿法的初始值为1。

三、实验仪器设备与材料

主流微型计算机

四、实验原理

二分法通过将含根区间逐步二分，从而将根的区间缩小到容许误差范围。牛顿通过迭代的方法逐步趋进于精确解，该两种方法的公式如下：

1、高斯-塞德尔迭代法

1）判断线性方程组是否主对角占优

n i a a

ii n

i

j j ij

,,2,1,1

=≤∑≠=

2）直接分离xi ，即

n i a x b d x ii n

j j ij i i ,,2,1,/)(1

=-=∑=

建立高斯-塞德尔迭代格式为：

n i a x a

x a d x ii n

i j k j

ij

i j k j

ij i k i

,,2,1,/)(1

)

(1

1

)

1()

1( =-

-=∑∑+=-=++

3）取初值迭代求解至所要求的精度为止。 2、牛顿法

)

()

(1k k k k x f x f x x '-

=+

五、实验步骤

1、理解并掌握二分法与牛顿法的公式；

2、画出二分法与牛顿法的流程图

3、使用C 语言编写出相应的程序并调试验证通过

六、实验报告要求

1、统一使用《武汉科技大学实验报告》本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。

2、源程序需打印后粘贴在实验报告册内；

3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。

七、实验注意事项

对于二分法应注意二分后如何判断根的区间，对于牛顿法注意如何确定迭代过程的结束

八、思考题

若使用牛顿法是发散的，如何对牛顿法进行改进以保证其收敛性。

前三个实验的程序代码（C/C++）和运行结果截图

Gauss 全选主元解方程组的源程序及运行结果

#include #include #include

using namespace std;

class Matrix {

public:

Matrix(); ~Matrix();

void SetMatrix(const int n,const double esp1);//构造线性方程组相应的矩阵，n 为方程的未知数数目，esp1为要求的精度 void Max(const int r);//全选主元

void ChangeRC(const int r);//根据主元变换矩阵的行或列

void Eliminate(const int r);//处理消元工作

void Result()const;//计算方程的解

void Calculate();

int GetRank()const;//返回矩阵的行数

double GetX(const i)const;//确定方程组的第i个解(1<=i<=N)

private:

//指针a和b分别用于存储方程组的未知数系数和方程"="右边的常数，esp存//储精度

double *a,*b,esp;

//指针flag用于记录方程组解的顺序

int *flag;

//以下的结构体用于在全选主元中记录最大主元的位置

struct coordinate

{

int row,column;

}location;

int N;//方程组未知数的数目

};

int main()

{

int n;

double esp1;

Matrix matrix;

do{

cout<<"请输入方阵的阶数：";

cin>>n;

if(n<0) n=0;//如何控制非法字符的输入？？？？？？？？？？

}while(n==0);

do

{

cout<<"请输入计算精度：";

cin>>esp1;

if(esp1<0) esp1=0;//输入不合法的精度就把精度置0

}while(esp1==0);

cout<<"输入线性方程组的增广矩阵:\n";

matrix.SetMatrix(n,esp1);//设置矩阵内的数据

matrix.Calculate();//计算方程组的解

//输出方程组的解

cout<<"\n\n方程组的解如下:\n";

for(int i=1;i<=matrix.GetRank();++i)

cout<<"X["<

return 0;

}

Matrix::Matrix()

{ //将Matrix类的数据成员初始化

a=NULL;

b=NULL;

flag=NULL;

location.row=0;

location.column=0;

esp=0;

N=0;

}

Matrix::~Matrix()

{ //释放指针a、b和flag指向的内存空间delete[]a;

delete[]b;

delete[]flag;

}

void Matrix::SetMatrix(const int n,const double esp1) {

N=n;

esp=esp1;

a=new double[N*N];

b=new double[N];

flag=new int[N];

//判断是否成功分配存储区

if(a==NULL||b==NULL||flag==NULL)

{

cout<<"分配存储区失败！\n";

exit(EXIT_FAILURE);

}

//读取线性方程组的增广矩阵

for(int i=0;i

{

for(int j=0;j>*(a+i*N+j);

cin>>*(b+i);

}

//flag中存储的值对应相应的x值，当方程的解由于列变换交换后，flag中//的值也相应交换，最后用于恢复解的顺序

for( i=0;i

}

void Matrix::Max(const int r)

{

double max=0;

for(int i=r;i

for(int j=r;j

{

if(max

{

max=fabs(*(a+i*N+j));

//设定最大主元的行、列

location.row=i;

location.column=j;

}

}

//最大主元小于输入的精度时，认为方程组无解，退出程序

if(max<=esp)

{

cout<<"方程组无解！\n";

exit(EXIT_FAILURE);

}

}

void Matrix::ChangeRC(const int r)

{

double temp;

//如果最大主元所在的行不在当前行，则进行行变换

if(location.row!=r)

{

for(int i=r;i

{

temp=*(a+r*N+i);

*(a+r*N+i)=*(a+location.row*N+i);

*(a+location.row*N+i)=temp;

}

temp=b[r];

b[r]=b[location.row];

b[location.row]=temp;

}

//若最大主元所在的列不在当前的r列，则进行列变换

if(location.column!=r)

{

for(int i=r;i

{

temp=*(a+i*N+r);

*(a+i*N+r)=*(a+i*N+location.column);

*(a+i*N+location.column)=temp;

}

//交换flag中的元素来标记方程解的位置变化

int temp1;

temp1=*(flag+r);

*(flag+r)=*(flag+location.column);

*(flag+location.column)=temp1;

}

}

void Matrix::Eliminate(const int r)

{

if(fabs(*(a+N*r+r))<=esp)

{

cout<<"方程组无解！\n";

exit(EXIT_FAILURE);

}

for(int i=r+1;i

{

for(int j=r+1;j

(*(a+i*N+j))-=(*(a+i*N+r))*(*(a+r*N+j))/(*(a+r*N+r));

(*(b+i))-=(*(b+r))*(*(a+i*N+r))/(*(a+r*N+r));

}

}

void Matrix::Result()const

{

if(fabs(*(a+N*(N-1)+N-1))<=esp)

{

cout<<"方程组无解！\n";

exit(EXIT_FAILURE);

}

double temp;

*(b+N-1)/=(*(a+N*(N-1)+N-1)); //求出X[N-1]

//依次求出X[i](i=N-2,N-3 (1)

for(int i=N-2;i>=0;--i)

{

temp=0;

for(int j=i+1;j

*(b+i)=(*(b+i)-temp)/(*(a+i*N+i));

}

//根据flag中的数据用冒泡排序法恢复方程组解的次序

for(i=0;i

for(int j=0;j

if(*(flag+j)>*(flag+j+1))

{

int temp1;

//交换解的顺序

temp=*(b+j);

*(b+j)=*(b+j+1);

*(b+j+1)=temp;

//交换用于标记的元素的顺序

temp1=*(flag+j);

*(flag+j)=*(flag+j+1);

*(flag+j+1)=temp1;

}

}

void Matrix::Calculate()

{ //根据矩阵行数重复进行寻找最大主元、变换行或列、消元for(int i=0;i

{

Max(i);

ChangeRC(i);

Eliminate(i);

}

Result();

}

int Matrix::GetRank()const

{

return N;//返回矩阵的行数

}

double Matrix::GetX(const int i)const

{

return *(b+i-1);

}

运行结果

追赶法求解方程组的算法：

1．输入方程组的维数n，将主对角元素b(i)(i=0:n-1)，，主对角元素左边的元素a(i)(i=0:n-2)，主对角元素右边的元素c(i)(i=0:n-2)，右端项的元素f(i)(i=0:n-1) 2．对方程组的系数矩阵作Crout分解, α(0)=b(0),对于i=0:n-2, c(i):=β(i):= c(i)/b(i), b(i+1):=α(i+1):= b(i+1)-a(i)* β(i)

3.解方程组Ly=f

b(0):=y(0):=[f(0)/ α(0)]:=[f(0)/b(0)]

对于i=1：n-1，b(i):=y(i):=[f(i)-a(i-1)*y(i-1)]/b(i)

4..解方程组Ux=y

a(n-1):=x(n-1):=b(n-1)

对于i=n-2:0,a(i)=x(i):=b(i)-c(i)*a(i+1);

5．输出方程组的解a(0:n-1)

用追赶法求解方程组的源程序及运行结果

#include

#include

using namespace std;

class MatrixThr

{

public:

MatrixThr();

~MatrixThr();

void SetMatrixThr(const int n);//设置三对角矩阵的数据

void Result();//计算三对角矩阵的解

double GetX(const int i)const;//取得第i个解,i从1开始

int GetN() const;//返回未知数的数目

private:

int N;//N为未知数的数目

//b为矩阵主对角线的元素首地址，a为主对角线左边一斜条元素的首地址，//c为主对角线右边一斜条元素首地址，f为方程组的常数首地址

double *a,*b,*c,*f;

};

int main()

{

MatrixThr matrix;

int n;

do{

cout<<"输入三对角方程组的变量的数目N:";

cin>>n;

}while(n<3);

cout<<"请依次输入三对角方程组每行的数据(0元素除外):\n";

matrix.SetMatrixThr(n); //计算方程组的解

matrix.Result();//输出方程组的解

cout<<"方程的解如下:\n";

for(int i=1;i<=matrix.GetN();++i)

cout<<"X["<

return 0;

}

MatrixThr::MatrixThr()

{ //初始化相关数据

N=0;

a=NULL;

b=NULL;

c=NULL;

f=NULL;

}

MatrixThr::~MatrixThr()

{ //释放分配的内存空间

delete []a;

delete []b;

delete []c;

delete []f;

}

void MatrixThr::SetMatrixThr(const int n)

{ //根据输入的未知数个数设置矩阵的数据

N=n;

a=new double[N];

b=new double[N];

c=new double[N];

f=new double[N];

//若内存分配失败，退出程序

if(a==NULL||b==NULL||c==NULL||f==NULL)

{

cout<<"初始化分配存储空间失败！\n";

exit(EXIT_FAILURE);

}

//依次输入三对角矩阵每行的元素

cin>>*b>>*c>>*f;

for(int i=1;i>*(a+i-1)>>*(b+i)>>*(c+i)>>*(f+i);

cin>>*(a+i-1)>>*(b+i)>>*(f+i);

}

void MatrixThr::Result()

{ //对系数矩阵A作Crout分解

for(int i=0;i

{ //将U中的α存于指针b中，L中的β存于指针c中*(c+i)/=(*(b+i));

*(b+i+1)-=(*(a+i))*(*(c+i));

}