398 lines
21 KiB
Fortran
398 lines
21 KiB
Fortran
! ==============================================================================
|
||
! 通过 gfortran ./test.f90 -o ./run && ./run 运行
|
||
! 程序名：
|
||
! 第11章习题
|
||
! 目的:
|
||
!
|
||
! 修订记录:
|
||
! 日期 编程者 改动描述
|
||
! =================== ============= =====================================
|
||
! 2021-05-13 14:32:35 Sola 习题11-5 判断格式是否正确
|
||
! 2021-05-13 15:42:31 Sola 习题11-6 函数的导数
|
||
! 2021-05-13 17:17:22 Sola 习题11-7 经时计算，判断单精度和双精度时间
|
||
! 2021-05-14 19:19:29 Sola 习题11-8 跳过，习题11-9 复数计算
|
||
! 2021-05-14 19:56:11 Sola 习题11-10 复数的振幅和相位
|
||
! 2021-05-14 20:09:50 Sola 习题11-11 欧拉公式
|
||
! 程序结构:
|
||
!
|
||
! ==============================================================================
|
||
! 模块：
|
||
module Chapter11
|
||
implicit none
|
||
! 数据字典
|
||
! 声明常数
|
||
REAL, PARAMETER :: PI=3.14159265 ! PI值
|
||
REAL, PARAMETER :: e=2.718281828459 ! 自然对数
|
||
INTEGER, PARAMETER :: arrayLength=20 ! 数组基准长度
|
||
! 声明变量
|
||
! 创建显式接口
|
||
contains
|
||
! subroutine SubName(varName1,varName2)
|
||
! implicit none
|
||
! ! 数据字典
|
||
! end subroutine SubName
|
||
end module Chapter11
|
||
! ==============================================================================
|
||
! 主程序：
|
||
program ProName
|
||
implicit none
|
||
! 数据字典
|
||
! 声明常量
|
||
! 声明变量
|
||
! 变量初始化
|
||
! 数据输入
|
||
! 运算过程
|
||
! 结果输出
|
||
! call Exercises11_5
|
||
! call Exercises11_6
|
||
! call Exercises11_7
|
||
! call Exercises11_9
|
||
! call Exercises11_10
|
||
call Exercises11_11
|
||
end program ProName
|
||
! ==============================================================================
|
||
! 子程序
|
||
! subroutine SubName(varName1,varName2)
|
||
! use MouName
|
||
! implicit none
|
||
! Type, intent(inout) :: varName
|
||
! end subroutine SubName
|
||
! ==============================================================================
|
||
! 函数
|
||
! function FunName(varName1,varName2)
|
||
! use MouName
|
||
! implicit none
|
||
! end function FunName
|
||
! ==============================================================================
|
||
! 习题11-5
|
||
subroutine Exercises11_5
|
||
implicit none
|
||
call Suba
|
||
! call Subb
|
||
contains
|
||
subroutine Suba
|
||
integer, parameter :: sgl = kind(0.0)
|
||
integer, parameter :: dbl = kind(0.0D0)
|
||
real(kind=sgl) :: a
|
||
real(kind=dbl) :: b
|
||
integer :: i
|
||
integer :: errorLevel
|
||
open(unit=1, status='scratch', iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
write(*, *) sgl, dbl
|
||
do i = 1, 45
|
||
write(*, *) i, selected_real_kind(r=i)
|
||
end do
|
||
write(1, '(A)') &
|
||
&"111111111111111111111111111111111111111111111", &
|
||
&"222222222222222222222222222222222222222222222"
|
||
rewind(unit=1, iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
read(1, '(F16.2)') a, b
|
||
write(*, *) a, b
|
||
end subroutine Suba
|
||
! subroutine Subb
|
||
! end subroutine Subb
|
||
end subroutine Exercises11_5
|
||
! 习题11-6 函数的导数
|
||
subroutine Exercises11_6
|
||
implicit none
|
||
integer, parameter :: realType = selected_real_kind(p=13) ! 确认实数类型
|
||
real(realType) :: x0, dx ! 待测点, 步长
|
||
real(realType) :: derivationX0 ! 导数值
|
||
real(realType), external :: Fx ! 外部函数
|
||
x0 = 0._realType ! 变量初始化
|
||
dx = 0.01_realType
|
||
call Derivation(Fx, x0, dx, derivationX0) ! 计算导数
|
||
write(*, *) '函数在0处的值为', Fx(0._realType), '在0处的导数值为', derivationX0 ! 输出结果
|
||
contains
|
||
subroutine Derivation(inputFunction, x0, dx, derivationX0)
|
||
integer, parameter :: realType = selected_real_kind(p=13) ! 确认实数类型
|
||
real(realType), intent(in) :: x0, dx ! 输入点位及步长
|
||
real(realType), intent(out) :: derivationX0 ! 输出的导数值
|
||
real(realType), external :: inputFunction ! 外部函数
|
||
derivationX0 = (inputFunction(x0 + dx) - inputFunction(x0))/dx ! 计算导数值
|
||
end subroutine Derivation
|
||
end subroutine Exercises11_6
|
||
function Fx(x)
|
||
implicit none
|
||
integer, parameter :: realType = selected_real_kind(p=13) ! 确认实数类型
|
||
real(realType), intent(in) :: x ! 输入点位
|
||
real(realType) :: Fx ! 定义函数输出类型
|
||
Fx = 10._realType*sin(20._realType*x) ! 计算函数值
|
||
end function Fx
|
||
! 习题11-7 经时计算 不过话说回来，现在的编译器在这边好像还是有优化的，，，结果可能不太准
|
||
subroutine Exercises11_7
|
||
implicit none
|
||
integer, parameter :: dbl = selected_real_kind(p=13)
|
||
integer, parameter :: sgl = selected_real_kind(p=1)
|
||
real(dbl), dimension(10, 10) :: matrix
|
||
real(dbl), dimension(10) :: arrayX
|
||
real(sgl), dimension(10, 10) :: matrix1
|
||
real(sgl), dimension(10) :: arrayX1
|
||
integer, dimension(8) :: timeNow
|
||
real :: timePast
|
||
integer, dimension(8) :: timeOld
|
||
integer :: errorLevel
|
||
integer :: i
|
||
open(unit=1, status='scratch', iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
write(1, '(A)') &
|
||
&" -2., 5., 1., 3., 4. -1., -2., -1., -5., -2.", &
|
||
&" 6., 4., -1., 6., -4., -5., 3., -1., 4., 3.", &
|
||
&" -6., -5., -2., -2., -3., 6., 4., 2., -6., 4.", &
|
||
&" 2., 4., 4., 4., 5., -4., 0., 0., -4., 6.", &
|
||
&" -4., -1., 3., -3., -4., -4., -4., 4., 3., -3.", &
|
||
&" 4., 3., 5., 1., 1., 1., 0., 3., 3., 6.", &
|
||
&" 1., 2., -2., 0., 3., -5., 5., 0., 1., -4.", &
|
||
&" -3., -4., 2., -1., -2., 5., -1., -1., -4., 1.", &
|
||
&" 5., 5., -2., -5., 1., -4., -1., 0., -2., -3.", &
|
||
&" -5., -2., -5., 2., -1., 3., -1., 1., -4., 4."
|
||
rewind(unit=1, iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
arrayX = (/ -5., -6., -7., 0., 5., -8., 1., -4., -7., 6./)
|
||
read(1, *) (matrix(i, :), i = 1, 10)
|
||
rewind(unit=1, iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
arrayX1 = (/ -5., -6., -7., 0., 5., -8., 1., -4., -7., 6./)
|
||
read(1, *) (matrix1(i, :), i = 1, 10)
|
||
call set_timer
|
||
do i = 1, 1000000
|
||
call GAEli1(matrix1, arrayX1, errorLevel)
|
||
! arrayX1 = (/ -5., -6., -7., 0., 5., -8., 1., -4., -7., 6./)
|
||
end do
|
||
call elapsed_time(timePast)
|
||
write(*, *) '使用单精度计算消耗时间', timePast, 's'
|
||
call set_timer
|
||
do i = 1, 1000000
|
||
call GAEli2(matrix, arrayX, errorLevel)
|
||
! arrayX = (/ -5., -6., -7., 0., 5., -8., 1., -4., -7., 6./)
|
||
end do
|
||
call elapsed_time(timePast)
|
||
write(*, *) '使用双精度计算消耗时间', timePast, 's'
|
||
contains
|
||
! 解方程子程序
|
||
! 经时子程序
|
||
subroutine set_timer ! 创建子程序1
|
||
call DATE_AND_TIME(values=timeNow) ! 调用当前时间子程序
|
||
end subroutine set_timer ! 结束子程序1
|
||
subroutine elapsed_time(timePast) ! 创建子程序2
|
||
real, intent(out) :: timePast ! 定义输出变量
|
||
timeOld = timeNow ! 传递值
|
||
call set_timer ! 调用子程序1
|
||
timePast = ((real(timeNow(3)-timeOld(3))*24 + real(timeNow(5)-timeOld(5)))&
|
||
&*60 + real(timeNow(6)-timeOld(6)))*60 + real(timeNow(7)-timeOld(7)) + &
|
||
&real(timeNow(8)-timeOld(8))/1000 ! 计算经历时间(秒)
|
||
end subroutine elapsed_time ! 结束子程序2
|
||
! 高斯-亚当消元法，不破坏输入矩阵
|
||
subroutine GAEli1(matrixInput1, arrayX, errorLevel)
|
||
integer, parameter :: varKind = selected_real_kind(p=2)
|
||
real(varKind), dimension(:, :), intent(inout) :: matrixInput1 ! 输入矩阵(输入)
|
||
real(varKind), dimension(size(matrixInput1, 1), size(matrixInput1, 2)) :: matrixInput
|
||
! 运算矩阵
|
||
real(varKind), dimension(:), intent(inout) :: arrayX ! 输入常数向量, 输出X解向量
|
||
integer, intent(out) :: errorLevel ! 错误值
|
||
real(varKind) :: temp ! 临时值, 用于数值交换
|
||
integer :: maxLocate ! 最大值所在行
|
||
integer :: i, n, m ! 局部变量: 循环参数
|
||
matrixInput = matrixInput1
|
||
m = size(matrixInput, 1)
|
||
do n = 1, m
|
||
maxLocate = n ! 初始化绝对值最大的位置
|
||
do i = n+1, m ! 执行最大支点技术，减少舍入误差
|
||
if ( abs(matrixInput(i, n)) > abs(matrixInput(maxLocate, n)) ) then
|
||
maxLocate = i
|
||
end if
|
||
end do
|
||
if ( abs(matrixInput(maxLocate, n)) < 1.0E-30 ) then
|
||
errorLevel = 1 ! 如果绝对值最大为0, 则矩阵不满秩, 无解
|
||
return ! 退出运算
|
||
end if
|
||
if ( maxLocate /= n ) then ! 如果绝对值最大值位置发生了变化, 则交换
|
||
do i = 1, m
|
||
temp = matrixInput(maxLocate, i)
|
||
matrixInput(maxLocate, i) = matrixInput(n, i)
|
||
matrixInput(n, i) = temp
|
||
end do
|
||
temp = arrayX(maxLocate)
|
||
arrayX(maxLocate) = arrayX(n)
|
||
arrayX(n) = temp
|
||
end if
|
||
do i = 1, m ! 消去其他行所有第n列的系数
|
||
if ( i /= n ) then
|
||
temp = -matrixInput(i, n)/matrixInput(n, n) ! 这里务必要设置一个临时值保存数值, 不然会被修改
|
||
matrixInput(i, :) = matrixInput(i, :) + matrixInput(n, :)*temp
|
||
arrayX(i) = arrayX(i) + arrayX(n)*temp
|
||
end if
|
||
end do
|
||
end do
|
||
do i = 1, m ! 产生X的解向量
|
||
arrayX(i) = arrayX(i)/matrixInput(i, i)
|
||
matrixInput(i, i) = 1.
|
||
end do
|
||
errorLevel = 0
|
||
end subroutine GAEli1
|
||
subroutine GAEli2(matrixInput1, arrayX, errorLevel)
|
||
integer, parameter :: varKind = selected_real_kind(p=13)
|
||
real(varKind), dimension(:, :), intent(inout) :: matrixInput1 ! 输入矩阵(输入)
|
||
real(varKind), dimension(size(matrixInput1, 1), size(matrixInput1, 2)) :: matrixInput
|
||
! 运算矩阵
|
||
real(varKind), dimension(:), intent(inout) :: arrayX ! 输入常数向量, 输出X解向量
|
||
integer, intent(out) :: errorLevel ! 错误值
|
||
real(varKind) :: temp ! 临时值, 用于数值交换
|
||
integer :: maxLocate ! 最大值所在行
|
||
integer :: i, n, m ! 局部变量: 循环参数
|
||
matrixInput = matrixInput1
|
||
m = size(matrixInput, 1)
|
||
do n = 1, m
|
||
maxLocate = n ! 初始化绝对值最大的位置
|
||
do i = n+1, m ! 执行最大支点技术，减少舍入误差
|
||
if ( abs(matrixInput(i, n)) > abs(matrixInput(maxLocate, n)) ) then
|
||
maxLocate = i
|
||
end if
|
||
end do
|
||
if ( abs(matrixInput(maxLocate, n)) < 1.0E-30 ) then
|
||
errorLevel = 1 ! 如果绝对值最大为0, 则矩阵不满秩, 无解
|
||
return ! 退出运算
|
||
end if
|
||
if ( maxLocate /= n ) then ! 如果绝对值最大值位置发生了变化, 则交换
|
||
do i = 1, m
|
||
temp = matrixInput(maxLocate, i)
|
||
matrixInput(maxLocate, i) = matrixInput(n, i)
|
||
matrixInput(n, i) = temp
|
||
end do
|
||
temp = arrayX(maxLocate)
|
||
arrayX(maxLocate) = arrayX(n)
|
||
arrayX(n) = temp
|
||
end if
|
||
do i = 1, m ! 消去其他行所有第n列的系数
|
||
if ( i /= n ) then
|
||
temp = -matrixInput(i, n)/matrixInput(n, n) ! 这里务必要设置一个临时值保存数值, 不然会被修改
|
||
matrixInput(i, :) = matrixInput(i, :) + matrixInput(n, :)*temp
|
||
arrayX(i) = arrayX(i) + arrayX(n)*temp
|
||
end if
|
||
end do
|
||
end do
|
||
do i = 1, m ! 产生X的解向量
|
||
arrayX(i) = arrayX(i)/matrixInput(i, i)
|
||
matrixInput(i, i) = 1.
|
||
end do
|
||
errorLevel = 0
|
||
end subroutine GAEli2
|
||
end subroutine Exercises11_7
|
||
! 习题11-9 复数计算
|
||
subroutine Exercises11_9
|
||
implicit none
|
||
integer, parameter :: dbl = selected_real_kind(p=3)
|
||
complex(dbl), dimension(3, 3) :: matrix
|
||
complex(dbl), dimension(3) :: arrayX
|
||
integer :: errorLevel
|
||
integer :: i
|
||
open(unit=1, status='scratch', iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
write(1, '(A)') &
|
||
&"( -2., 5.), ( 1., 3.), ( 4., -1.)", &
|
||
&"( 2., -1.), ( -5., -2.), ( -1., 6.)", &
|
||
&"( -1., 6.), ( -4., -5.), ( 3., -1.)", &
|
||
&"( 7., 5.), (-10., -8.), ( -3., -3.)"
|
||
rewind(unit=1, iostat=errorLevel)
|
||
if ( errorLevel /= 0 ) stop ""
|
||
read(1, *) (matrix(i, :), i = 1, 3)
|
||
read(1, *) arrayX
|
||
call GAEli(matrix, arrayX, errorLevel)
|
||
write(*, *) '计算结果为：'
|
||
write(*, *) 'x1 = ', arrayX(1)
|
||
write(*, *) 'x2 = ', arrayX(2)
|
||
write(*, *) 'x3 = ', arrayX(3)
|
||
contains
|
||
! 复数求解
|
||
subroutine GAEli(matrixInput1, arrayX, errorLevel)
|
||
integer, parameter :: varKind = selected_real_kind(p=3)
|
||
complex(varKind), dimension(:, :), intent(inout) :: matrixInput1 ! 输入矩阵(输入)
|
||
complex(varKind), dimension(size(matrixInput1, 1), size(matrixInput1, 2)) :: matrixInput
|
||
! 运算矩阵
|
||
complex(varKind), dimension(:), intent(inout) :: arrayX ! 输入常数向量, 输出X解向量
|
||
integer, intent(out) :: errorLevel ! 错误值
|
||
complex(varKind) :: temp ! 临时值, 用于数值交换
|
||
integer :: maxLocate ! 最大值所在行
|
||
integer :: i, n, m ! 局部变量: 循环参数
|
||
matrixInput = matrixInput1
|
||
m = size(matrixInput, 1)
|
||
do n = 1, m
|
||
maxLocate = n ! 初始化绝对值最大的位置
|
||
do i = n+1, m ! 执行最大支点技术，减少舍入误差
|
||
if ( cabs(matrixInput(i, n)) > cabs(matrixInput(maxLocate, n)) ) then
|
||
maxLocate = i
|
||
end if
|
||
end do
|
||
if ( cabs(matrixInput(maxLocate, n)) < 1.0E-30 ) then
|
||
errorLevel = 1 ! 如果绝对值最大为0, 则矩阵不满秩, 无解
|
||
return ! 退出运算
|
||
end if
|
||
if ( maxLocate /= n ) then ! 如果绝对值最大值位置发生了变化, 则交换
|
||
do i = 1, m
|
||
temp = matrixInput(maxLocate, i)
|
||
matrixInput(maxLocate, i) = matrixInput(n, i)
|
||
matrixInput(n, i) = temp
|
||
end do
|
||
temp = arrayX(maxLocate)
|
||
arrayX(maxLocate) = arrayX(n)
|
||
arrayX(n) = temp
|
||
end if
|
||
do i = 1, m ! 消去其他行所有第n列的系数
|
||
if ( i /= n ) then
|
||
temp = -matrixInput(i, n)/matrixInput(n, n) ! 这里务必要设置一个临时值保存数值, 不然会被修改
|
||
matrixInput(i, :) = matrixInput(i, :) + matrixInput(n, :)*temp
|
||
arrayX(i) = arrayX(i) + arrayX(n)*temp
|
||
end if
|
||
end do
|
||
end do
|
||
do i = 1, m ! 产生X的解向量
|
||
arrayX(i) = arrayX(i)/matrixInput(i, i)
|
||
matrixInput(i, i) = 1.
|
||
end do
|
||
errorLevel = 0
|
||
end subroutine GAEli
|
||
end subroutine Exercises11_9
|
||
! 习题11-10 复数的振幅和相位
|
||
subroutine Exercises11_10
|
||
implicit none
|
||
integer, parameter :: sgl=selected_real_kind(p=2) ! 定义单精度
|
||
complex(sgl) :: var1 ! 定义复数
|
||
real(sgl) :: amp, theta ! 定义振幅和相位
|
||
real(sgl), parameter :: PI=3.14159265 ! 圆周率
|
||
call InputComplex(amp, theta) ! 调用子程序, 获取输入复数的振幅和相位
|
||
write(*, *) '振幅为：', amp, '， 相位为：', theta, '°' ! 输出结果
|
||
contains
|
||
subroutine InputComplex(amp, theta)
|
||
integer, parameter :: sgl=selected_real_kind(p=2) ! 定义单精度
|
||
real(sgl), intent(out) :: amp, theta ! 定义输出结果
|
||
complex(sgl) :: var1 ! 定义输入复数
|
||
write(*, *) '请输入一个复数：' ! 提示信息
|
||
read(*, *) var1 ! 读取输入复数
|
||
amp = cabs(var1) ! 获取振幅
|
||
theta = atan(aimag(var1)/real(var1))*360./(2.*PI) ! 获取相位(角度)
|
||
end subroutine InputComplex
|
||
end subroutine Exercises11_10
|
||
! 习题11-11 欧拉公式
|
||
subroutine Exercises11_11
|
||
implicit none
|
||
integer, parameter :: sgl=selected_real_kind(p=2) ! 定义单精度
|
||
real(sgl) :: theta ! 定义角度(弧度)
|
||
complex(sgl) :: var1 ! 定义复数
|
||
real(sgl) :: PI=3.14159265 ! 圆周率
|
||
integer :: i ! 循环参数
|
||
do i = 0, 2 ! 循环theta = 0, pi/2, pi
|
||
theta = i*PI/2. ! 计算theta
|
||
write(*, '("theta = ", F6.2,", e^theta_i = " 2("(", ES9.2, ",", ES9.2, ") "))') &
|
||
&theta, cexp(cmplx(0., theta, sgl)), EulerFormula(theta)! 输出结果
|
||
end do
|
||
contains
|
||
! 欧拉公式
|
||
function EulerFormula(theta)
|
||
integer, parameter :: sgl=selected_real_kind(p=2) ! 定义单精度
|
||
real(sgl), intent(in) :: theta ! 定义输入角度(弧度)
|
||
complex(sgl) :: EulerFormula ! 定义函数返回类型
|
||
EulerFormula = cmplx(cos(theta), sin(theta), sgl) ! 计算返回值
|
||
end function EulerFormula
|
||
end subroutine Exercises11_11 |