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;ifor(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;ifor(intj=0;jif(*(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