fortran数值计算基础.docx

fortran数值计算基础.docx

《fortran数值计算基础.docx》由会员分享，可在线阅读，更多相关《fortran数值计算基础.docx（59页珍藏版）》请在冰豆网上搜索。

fortran数值计算基础

数值计算根底

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

实验二插值方法11

实验三数值积分5

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

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

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

一、实验目的

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

二、实验内容

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

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

2.5

2.3

-5.1〕

■xj

3.71

5.3

9.6

1.5

X2

=3.8

8.1

1.7

—4.3_

[X3_

虬5_

2、用追赶法求解方程组

_-20

0

0

01

一X」

一-10〕

1

-2

0

0

0

X2

0

0

1

-2

0

0

X3

—

0

0

0

1

-2

0

X4

0

■.0

0

0

1

-2」

*一

1

0一

三、实验仪器设备与材料

主流微型电脑

四、实验原理

1、Guass列选主元消去法

对于AX=B

~~~

1）、消元过程：

将（A|B）进行变换为了（A|B）,其中A是上二角矩阵.即:

■A

a〔1a〔2a〔n

b1

1a〔2a〔n

b1

a21a22a2n

b2

01…a2n

b2

aaa

a

T

+盲a

a

・—■

■

——■

■

■■■

■

■■■

■

an1an2…ann

bnJ

<00…ann

bn/

k从1到n-1

a、列选主元

选取第k列中绝对值最大元素maxlaikl作为了主元.

b、换行

akj：

二aj,j=k1,,n

c、回一化

akj/akk=akj,j=k1,,nbk/akk=bk

d、消元

aj—aikakj=a0,i=k1,,n;j=k1,,nbifbk=bi,i=k1,,n

~一~一

2）、回代过程:

由（A|B）解出xn,xn」,’,x1.

bn/ann='

Xn

n

bk-'akjXj=Xk,k=n-1,,2,1

j=k1

2、追赶法

线性方程组为了：

角

G

ZX1'

p1\

b2

a2

C2

X2

f2

b3

a3

+

C3

+

+

X3

■■

—

f3

♦.

+

bnJ

+

an」

Cn^L

■■

Xna

■■

bn

anj

5」

做LU分解为了:

『％

L«2

+

分解公式:

%=ai（i=2,3,…,n）

皿=bi,M=biW—（i=2,3,…,n）

穴=色（i=1,2,…,n—1）:

i

那么

"Ly=f

Ax=fnLUx=fn』Ux=y

回代公式：

―fi

yi=%

f

y'''（i=2,3,,n）

=yi—EiX.（i=n—1,n—2,…,1）

五、实验步骤

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

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

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

六、实验报告要求

1、统一使用〈〈武汉科技大学实验报告>本书写,实验报告的内容要求有：

实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个局部.

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

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

七、实验考前须知

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

八、思考题

假设使用全主元消去法,在编程中应如何记录保存对于未知数的调换.

实验二插值方法

一、实验目的

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

二、实验内容

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

以下函数表

Xi

0.56i60

0.56280

0.5640i

0.5652i

y

0.8274i

0.82659

0.82577

0.82495

求x=0.5635时的函数值.

三、实验仪器设备与材料

主流微型电脑

四、实验原理

n个插值节点的函数值,那么可由拉格郎日插值公式与牛顿插值公式构造出插值多项式,从而由该插值多项式求出所要求点的函数值.拉格郎日插值公式与牛顿插值公式如下：

1、Lagrange插值公式

n

Ln（x）=l°（x）y0li（x）yi...Tn（x）yn='y」k（x）

k=0

nx-xj

=H

jTXk_Xjj=k

lk（x）=（x-x0）（x-xi）（X-Xk^Xfi）（XE（Xk—X°）（Xk—Xi）（Xk—X“）（Xk—Xki）（Xk—Xn）

2、Newton插值公式

Nn（x）=f（Xo）f[X0,Xi]（X-X0）f[X0,Xi,X2]（X-Xo）（X-Xi）f[X0,Xi,Xn]（X-X°）（X-Xi广（X-Xna）五、实验步骤

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

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

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

六、实验报告要求

i、统一使用〈〈武汉科技大学实验报告>本书写,实验报告的内容要求有：

实验目的、

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

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

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

七、实验考前须知

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

八、思考题

比拟Lagrange插值法与Newton插值法的异同.

实验三数值积分

一、实验目的

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

、实验内容

分别写出变步长梯形法与Romberge法计算定积分的算法,编写程序上机调试出结果,

要求所编程序适合于任何类型的定积分,即能解决这一类问题,而不是某一个问题.实验中

以以下数据验证程序的正确性.

炎尝dx,&<0.00001o0x

三、实验仪器设备与材料

主流微型电脑

四、实验原理

通过变步长梯形法与龙贝格法,我们只要知道n个求积节点的函数值,那么可由相应

的公式求出该函数的积分值,从而不需要求该函数的原函数.变步长梯形法与龙贝格法公式如下：

1、变步长梯形法

n1h

Tn八h[f（Xi）f（Xi）]

1=02

hn4

—[f（a）2、5）f（b）]

2i4

n』

2i=0

f（xi1/2）

用丁2n—Tn|专&来控制精度

2、龙贝格法

丁2"=上+壬f（xf/2）

22i*

Sn

Cn

Rn

=4T-1T

cT2ncTn

33

16c1c

S2^—Sn

1515

=^C2n—Cn

6363

R2n-Rn

＜耳来控制精度

五、实验步骤

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

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

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

六、实验报告要求

1、统一使用〈〈武汉科技大学实验报告>本书写,实验报告的内容要求有：

实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个局部.

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

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

七、实验考前须知

isinx

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

的定义.

八、思考题

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

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

一、实验目的

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

二、实验内容

分别写出改良欧拉法与四阶龙格-库塔求解的算法,编写程序上机调试出结果,要求所

编程序适合于任何一阶常微分方程的数值解问题,即能解决这一类问题,而不是某一个问题.

实验中以以下数据验证程序的正确性.

2

求厂——约步长h=0.25.

、y（0）=2（0

三、实验仪器设备与材料

主流微型电脑

四、实验原理

常微分方程的数值解主要采用“步进式〞,即求解过程顺着节点排列次序一步一步向前

推进,在单步法中改良欧拉法和四阶龙格-库塔法公式如下：

1、改良欧拉法

yn1=ynhf（Xn,yn）

h,、一

y1=yn2[f（Xn,Yn尸f（Xn1,yn1）]

2、四阶龙格-库塔法

ki

yn1=yn*k12k22k‘kQ

=f（Xn,yn）

=f（Xn^,yn-hk1）

22

k3

hh..

=f（Xn,ynk2）22

=f（Xnh,ynhk3）

五、实验步骤

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

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

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

六、实验报告要求

1、统一使用〈〈武汉科技大学实验报告>本书写,实验报告的内容要求有：

实验目的、

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

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

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

七、实验考前须知

2

"=一'"的精确解为了y=2/（1十x2）,通过调整步长,打量结果的

、y（0）=2（0^x壬5）

精度的改变

八、思考题

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

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

一、实验目的

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

二、实验内容

分别写出高斯-塞德尔迭代法与牛顿迭代法的算法,编写程序上机调试出结果,要求所

编程序适合于任何一个方程的求根,即能解决这一类问题,而不是某一个问题.实验中以下

列数据验证程序的正确性.

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

一7

21

-2]"

一4]

9

153

—2x2

7

-2

-211

5x3

--1

_1

32

13一为了一

一.一

2、用牛顿迭代法求方程x‘—x—1=0的近似根,&<0.00001,牛顿法的初始值为了1.

三、实验仪器设备与材料

主流微型电脑

四、实验原理

二分法通过将含根区间逐步二分,从而将根的区间缩小到容许误差范围.牛顿通过迭代

的方法逐步趋进于精确解,该两种方法的公式如下：

1、高斯-塞德尔迭代法

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

n

E

j任j#

aij

<

aii

i=1,2,…,n

2）直接别离

xi,即

xi

=（di

n

-Z

bijxj）/a,i=1,2,…,n

j1

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

i-1

（k1）（k1）

xi-（di―aijxj

j4

n

-'aijXj〞））/aii,i=1,2,,n

j壬1

3）取初值迭代求解至所要求的精度为了止.

2、牛顿法

Xs=Xk

f（Xk）

f（Xk）

五、实验步骤

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

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

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

六、实验报告要求

1、统一使用〈〈武汉科技大学实验报告>本书写,实验报告的内容要求有：

实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个局部.

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

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

七、实验考前须知

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

八、思考题

假设使用牛顿法是发散的,如何对牛顿法进行改良以保证其收敛性.

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

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

#include

#include

#include

usingnamespacestd;

classMatrix

{

public:

Matrix（）;

~Matrix（）;

voidSetMatrix（constintn,constdoubleesp1）;/构造线性方程组相应的矩阵,

为了方程的未知数数I,esp1为了要求的精度

voidMax（constintr）;//全选主元

voidChangeRC（constintr）;//恨据主元变换矩阵的行或歹U

voidEliminate（constintr）;〃处理消元工作

voidResult（）const;//计算方程的解

voidCalculate（）;

intGetRank（）const;//返回矩阵的行数

doubleGetX（consti）const;溯定方程组的第i个解（1<=i<=N）

private:

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

〃储精度

double*a,*b,esp;

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

int*flag;

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

structcoordinate

{

introw,column;

}location;

intN;//方程组未知数的数目

};

intmain（）

{

intn;

doubleesp1;

Matrixmatrix;

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（inti=1;i<=matrix.GetRank（）;++i）

cout<<"X["<

"<

return0;

}

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;

}voidMatrix:

:

SetMatrix（constintn,constdoubleesp1）{

N=n;

esp=esp1;

a=newdouble[N*N];

b=newdouble[N];

flag=newint[N];

〃判断是否成功安排存储区

if（a==NULL||b==NULL||flag==NULL）

{

cout<<"安排存储区失败！

\n";

exit（EXIT_FAILURE）;

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

for（inti=0;i

（

for（intj=0;j>*（a+i*N+j）;

cin>>*（b+i）;

}

//flag中存储的值对应相应的x值,当方程的解由于列变换交换后,flag中

〃的值也相应交换,最后用于恢复解的顺序

for（i=0;i

}

voidMatrix:

:

Max（constintr）

（

doublemax=0;

for（inti=r;i

for（intj=r;j

（

if（max

（

max=fabs（*（a+i*N+j））;

〃设定最大主元的行、列

location.row=i;

location.column=j;

}

}

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

if（max<=esp）

（

cout<<"方程组无解！

\n〞;

exit（EXIT_FAILURE）;

}

}voidMatrix:

:

ChangeRC（constintr）

（

doubletemp;

〃如果最大主元所在的行不在当前行,那么进行行变换

if（location.row!

=r）

（

for（inti=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（inti=r;i

{

temp=*（a+i*N+r）;

*（a+i*N+r）=*（a+i*N+location.column）;

*（a+i*N+location.column）=temp;

}

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

inttemp1;

temp1=*（flag+r）;

*（flag+r）=*（flag+location.column）;

*（flag+location.column）=temp1;

}

}voidMatrix:

:

Eliminate（constintr）

{

if（fabs（*（a+N*r+r））<=esp）

{

cout<<"方程组无解！

\n";

exit（EXIT_FAILURE）;

}

for（inti=r+1;i

{

for（intj=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））;

}

}voidMatrix:

:

Result（）const

if（fabs（*（a+N*（N-1）+N-1））<=esp）

（

cout<<"方程组无解！

\n";

exit（EXIT_FAILURE）;

}

doubletemp;

*（b+N-1）/=（*（a+N*（N-1）+N-1））;//求出X[N-1]

〃依次求出X[i]（i=N-2,N-3•…,1）

for（inti=N-2;i>=0;--i）

（

temp=0;

for（intj=i+1;j

*（b+i）=（*（b+i）-temp）/（*（a+i*N+i））;

}

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

for（i=0;i

for（intj=0;j

if（*（flag+j）>*（flag+j+1））

（

inttemp1;

〃交换解的顺序

temp=*（b+j）;

*（b+j）=*（b+j+1）;

*（b+j+1）=temp;

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

temp1=*（flag+j）;

*（flag+j）=*（flag+j+1）;

*（flag+j+1）=temp1;

}

}voidMatrix:

:

Calculate（）

（〃根据矩阵行数重复进行寻找最大主元、变换行或列、消元

for（inti=0;i

（

Max（i）;

ChangeRC（i）;

Eliminate（i）;

}

Result（）;

intMatrix:

:

GetRank（）const

（

returnN;//返回矩阵的行数

}

doubleMatrix:

:

GetX（constinti）const

（

return*（b+i-1）;

}

运行结果

9.61.53.8

B»JL1.-4.35-.E-

方程组的解如下,

mi*.箱g河

-A,419325

I'l'cas云n亨yt：

.Qonii;inu.ie

半:

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

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分解,a（0）=b（0）,对丁i=0:

n-2,

c（i）:

=6（i）:

=c（i）/b（i）,b（i+1）:

=a（i+1）:

=b（i+1）-a（i）*6（i）

3.解方程组Ly=f

b（0）:

=y（0）:

=[f（0）/a（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

usingnamespacestd;

classMatri