利用中心差分格式数值求解导数讲课讲稿.docx
《利用中心差分格式数值求解导数讲课讲稿.docx》由会员分享，可在线阅读，更多相关《利用中心差分格式数值求解导数讲课讲稿.docx（13页珍藏版）》请在冰豆网上搜索。
利用中心差分格式数值求解导数讲课讲稿
利用中心差分格式数值求解导数
利用中心差分格式数值求解导数
一、问题描述
利用中心差分格式近似导数
，数值求解
步长分别取
二、格式离散
将x轴上[0,1]之间的线段按上述步长，等步长的离散为n个小段，包括端点，共n+1个网格节点，示意图如下：
线段上边的数字表示x轴上的坐标值，线段下边的数字表示节点编号，从0到n编号。
二阶导数中心差格式离散
整理为线性方程形式
其中，
为空间离散步长；i=1,2，……，n-1
包括边界条件的线性方程组如下：
改写成矩阵形式：
其中，
，
，
系数矩阵A中仅三对角线上的数值不全为0，其余位置上的数值全为0，是典型的对角占优的三对角矩阵，列向量f中，
，且
，作为边界条件。
追赶法求解线性方程组简述
对A做Crout分解，即
，
其中
为待定常数，由矩阵乘法可以得到下面的式子：
将对角占优三对角矩阵线性方程组
等价为如下两个方程组
，
求解对角占优三对角矩阵线性方程组的追赶法步骤：
①输入数据
②计算
③求解方程组
④求解方程组
⑤输出
计算流程图
三、程序中主要符号和数组意义
符号或数组
意义
A、B、C、D、filename
用于自动更改dat文件名的字符串变量
h
离散步长
n
离散网格数，共n+1个网格节点
p
辅助变量，暂时记录网格节点上的y值
数组x，y
离散节点的x，y坐标
子程序数组a,b,c
记录系数矩阵占优对角线上的值
子程序数组f
记录线性方程组常数向量
子程序数组s,r,t,g
追赶法求解线性方程组需要用到的中间辅助向量
四、计算结果与讨论
不同步长的数值结果函数曲线与精确解的对比
从对比图中可以看到，在所取的四个计算步长下，数值计算结果与精确解都符合得相当好
Step
Accuracy（Max-error）
0.05
8.152404966677018E-005
0.01
3.264604683472783E-006
0.001
3.817221055912867E-008
0.0001
9.862256566961491E-009
不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量，即上表中的Max-error来度量。
从上表中可以看出，随着网格节点的加密，Max-error的数量级在降低，即数值解的精度提高。
上述四个步长中，将线段离散成一万个网格时，数值结果的精度最高。
五、源程序
programmain
implicitnone
character（13）filename!
定义了五个字符串变量，用于按输出数据的需要自动更改dat文件的文件名
character（3）A
character（6）B
character（4）C
character（3）D
integer:
:
n,i
real（8）:
:
h,error,p
real（8）,allocatable:
:
y（:
）,x（:
）!
声明可变长度数组，x、y轴坐标定义为可变数组，数组长度按需要自动更改
write（*,*）"输入步长：
"
read（*,*）h!
读入空间步长
n=NINT（1.0/h）!
nint命令，取与浮点数最接近的整数
allocate（y（0:
n）,x（0:
n））!
指定可变数组的长度
A="po-"!
po代表problem
open（unit=11,file="help.txt"）!
打开一个txt文件，用于帮助更改dat文件的文件名
write（11,"（f6.4）"）h
rewind（11）!
将文件的读写位置移回到文件的最前面
read（11,"（A6）"）B
C=".dat"
filename=A//B//C
callsubsolution（y,n,h）!
调用追赶法求解线性方程组的子程序
open（unit=10,file=filename）
doi=0,n
x（i）=real（i）*h
enddo
write（10,*）'TITLE="X-YCURVE"'!
写入到dat文件中的一段字符，便于用tecplot软件后处理计算数据
write（10,*）'VARIABLES="X","Y"'
write（10,"（'ZONET=""Problem1-',f8.5,'"",I=',I6,',F=POINT'）"）h,n+1
doi=0,n
write（*,"（F10.8,10x,F10.8）"）x（i）,y（i）
write（10,"（F10.8,10x,F10.8）"）x（i）,y（i）!
将数值解数据写入dat文件
enddo
y（0）=0
y（n）=1
error=0.0
doi=1,n-1
p=y（i）
y（i）=（1+0.25*sin（2.0））*（i*h）-0.25*sin（2.0*i*h）!
求解节点精确值
error=max（error,abs（p-y（i）））!
寻找各个节点误差的最大值
enddo
write（10,"（'ZONET=""Problem1-',f8.5,'-exact"",I=',I6,',F=POINT'）"）h,n+1
doi=0,n
write（*,"（F10.8,10X,F10.8）"）x（i）,y（i）
write（10,"（F10.8,10X,F10.8）"）x（i）,y（i）!
将精确解数据写入dat文件
enddo
D="er-"!
er代表error，这里指精确值和数值计算值之间的差别
filename=D//B//C
open（unit=12,file=filename）
write（12,*）"max-error=",error!
将误差最大值写入dat文件
write（*,*）n,error
write（*,*）filename
stop
end!
主程序结束
subroutinesubsolution（y,n,h）!
子程序
implicitnone
integer:
:
n,i
real（8）:
:
h
real（8）:
:
y（0:
n）!
数组和变量的声明不能同时进行
real（8）,allocatable:
:
a（:
）,b（:
）,c（:
）,s（:
）,t（:
）,f（:
）,g（:
）!
声明可变长度数组
allocate（a（1:
n）,b（0:
n）,c（0:
n-1）,s（0:
n）,t（0:
n-1）,f（0:
n）,g（0:
n））!
指定可变数组的长度
a=1!
对数组a，b，c赋值
b=-2
c=1
a（n）=0
b（0）=1
b（n）=1
c（0）=0
f（0）=0!
对数组f赋值
f（n）=1
doi=1,n-1
f（i）=h*h*sin（2.0*real（i）*h）
enddo
s（0）=b（0）!
计算s，t
t（0）=c（0）/b（0）
doi=1,n-1
s（i）=b（i）-a（i）*t（i-1）
t（i）=c（i）/s（i）
enddo
s（n）=b（n）-a（n）*t（n-1）
g（0）=f（0）/b（0）!
追过程
doi=1,n
g（i）=（f（i）-a（i）*g（i-1））/s（i）
enddo
y（n）=g（n）!
赶过程
doi=n-1,1,-1
y（i）=g（i）-t（i）*y（i+1）
enddo
y（0）=0
y（n）=1
return
endsubroutinesubsolution
升级会员 