Sola/Fortran-95-2003-Program-3rd-Edition-Exercise
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update file date

Browse Source
This commit is contained in:
SolaSUSTech
2025-09-25 16:46:47 +08:00
parent 88fdb762e0
commit 1d9f18d8f3
15 changed files with 0 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

517
20210426-第8章习题.f90 Normal file
View File

@ -0,0 +1,517 @@
! ==============================================================================
! 通过 gfortran ./test.f90 -o ./run && ./run 运行
! 程序名：
! 习题
! 目的:
! 练习用
! 修订记录:
! 日期 编程者 改动描述
! =================== ============= =====================================
! 2021-04-26 16:53:03 Sola 测验8-1 7 测试数组调用
! 2021-04-30 18:03:39 Sola 习题8-7 输出数组中正负数和0的个数
! 2021-04-30 18:49:30 Sola 习题8-8 可分配数组输入数据定义并计算行列和
! 2021-04-30 19:26:59 Sola 习题8-10 没看懂啥意思, 跳过
! 2021-04-30 19:27:16 Sola 习题8-11 在赋值中修改数组范围，太长了，跳过
! 2021-04-30 19:30:55 Sola 习题8-12 DO循环和WHERE结构比较
! 2021-04-30 20:10:23 Sola 习题8-13 计算年平均温度
! 2021-04-30 20:46:31 Sola 习题8-14 矩阵乘法 8-15 8-16
! 2021-05-01 00:32:18 Sola 习题8-17 相对极大(鞍点)
! 2021-05-01 02:01:58 Sola 习题8-18 金属盘温度分布
! 程序结构:
!
! ==============================================================================
! 模块：
module Module_8
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 Module_8
! ==============================================================================
! 主程序：
program ProName
implicit none
! call Exercises8_1_7
! call Exercises8_7
! call Exercises8_8
! call Exercises8_9
! call Exercises8_10
! call Exercises8_11
! call Exercises8_12
! call Exercises8_13
! call Exercises8_14
! call Exercises8_17
call Exercises8_18
! 数据字典
! 声明常量
! 声明变量
! 变量初始化
! 数据输入
! 运算过程
! 结果输出
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
! ==============================================================================
! 测试8-1.7 使用数组从数组中调取元素
subroutine Exercises8_1_7
implicit none
integer, dimension(4) :: list1=(/1,4,2,2/) ! 定义数组1
integer, dimension(3) :: list2=(/1,2,3/) ! 定义数组2
integer, dimension(5, 5) :: array ! 定义数组3
integer :: i, j ! 循环参数
do i = 1, 5
do j = 1, 5
array(i, j) = i + 10*j ! 数组3赋值
end do
end do
write (*,'(1X, 4I4)') array(list1, list2) ! 测试通过数组调用数组
end subroutine Exercises8_1_7
! 习题8-7 统计数组中正负数和0的个数
subroutine Exercises8_7
implicit none
real, dimension(-50:50, 0:100) :: values ! 定义随机数组
integer :: i, j ! 循环参数
integer :: numPositive, numNegative, numZero ! 统计参数
! 初始化变量
numPositive = 0
numNegative = 0
numZero = 0
! FORALL方式随机赋值(好像不能这么用, 所以注释)
! forall (i=-50:50, j=0:100)
! call random_seed()
! call random_number(values(i, j))
! end forall
! 嵌套循环赋值
do i = -50, 50
do j = 0, 100
call random_seed() ! 根据日期, 时间获取随机数种子
call random_number(values(i, j)) ! 调用随机数子程序给数组赋值
end do
end do
values = values * 200. - 100. ! 调整数组范围
! 循环计算正负数和0的个数
do i = -50, 50
do j = 0, 100
if (values(i, j) > 0) then
numPositive = numPositive + 1
elseif (values(i, j) < 0) then
numNegative = numNegative + 1
else
numZero = numZero + 1
end if
end do
end do
! 输出结果
write(*,1) ' 测试组', numPositive, numNegative, numZero
write(*,1) ' 对照组', count(values > 0.), count(values < 0.), count(values == 0.)
1 format(1X, A/&
&1X, '数组中有正数 ', I5, '个'/&
&1X, '数组中有负数 ', I5, '个'/&
&1X, '数组中有0 ', I5, '个')
end subroutine Exercises8_7
! 习题8-8 可分配数组输入数据定义并计算行列和, 注意IO操作的返回值参数和数组分配不同
subroutine Exercises8_8
implicit none
real, dimension(:, :), allocatable :: array ! 定义可分配数组
integer :: errorLevel ! 错误码
integer :: a, b ! 数组范围
integer :: i ! 循环参数
character(len=20) :: inforError ! 错误信息
! 打开一个临时文件
open(unit=1, status='scratch', iostat=errorLevel, iomsg=inforError)
if (errorLevel /= 0) stop inforError ! 检测打开文件是否成功
! 输入参数 8-8
! write(1, *) 2, 4
! write(1, *) -24., -1121., 812.1, 11.1
! write(1, *) 35.6, 8.1E3, 135.23, -17.3
! 输入参数 8-9
write(1, *) 4, 5
write(1, *) 33., -12., 16., 0.5, -1.9
write(1, *) -6., -14., 3.5, 11., 2.1
write(1, *) 4.4, 1.1, -7.1, 9.3, -16.1
write(1, *) 0.3, 6.2, -9.9, -12., 6.8
! write(1, *) , , , ,
rewind(unit=1) ! 回到文件开头
read(1,*) a, b ! 读取数组范围
! 检查数组范围
if (a > 100) a = 100
if (b > 100) b = 100
allocate(array(a, b), stat=errorLevel, errmsg=inforError) ! 分配数组
if (errorLevel /= 0) stop inforError ! 检测数组分配是否成功
do i = 1, a
read(1, *) array(i, :) ! 数组赋值
end do
! 输出统计结果
1 format(1X, 'Sum of ', A, ' ', I3, ' = ', F10.3)
do i = 1, a
write(*, 1) 'row', i, sum(array(i, :))
end do
do i = 1, b
write(*, 1) 'col', i, sum(array(:, i))
end do
if (allocated(array)) deallocate(array, stat=errorLevel)
if (errorLevel /= 0) print *, "array: Deallocation request denied"
close(unit=1)
end subroutine Exercises8_8
! 习题8-10 没看懂啥意思, 跳过
! 习题8-11 在赋值中修改数组范围，太长了，跳过, 写了个简化版的，利用F03的特性实现
! 不过有一点需要注意，好像不可以用allocate语句多次定义，有点理解无缝是什么意思了
subroutine Exercises8_11
implicit none
real, dimension(:), allocatable :: array ! 定义可分配数组
integer :: i ! 循环变量
integer :: errorLevel ! 错误码
character(len=20) :: inforError ! 错误信息
real :: temp ! 临时变量
allocate(array(1), stat=errorLevel, errmsg=inforError) ! 分配数组
if (errorLevel /= 0) stop inforError ! 检测数组分配状态
call random_seed() ! 根据日期时间产生随机数种子
call random_number(array(1)) ! 用随机数给数组赋值
do i = 2, 50 ! 循环产生随机数, 并给数组赋值
call random_seed()
call random_number(temp)
array = (/array, temp/) ! 利用F03特性扩展数组
end do
write(*,'(1X, 10F6.3, " ")') array ! 10个一行输出数组内容
if (allocated(array)) deallocate(array, stat=errorLevel)
if (errorLevel /= 0) print *, "array: Deallocation request denied"
end subroutine Exercises8_11
! 习题8-12 DO循环和WHERE结构比较
subroutine Exercises8_12
implicit none
integer, parameter :: arrayLength=10 ! 设置数组基本长度
real, dimension(100*arrayLength, arrayLength, 3*arrayLength) :: arr, arrayTemp ! 定义数组
integer :: i, j, k ! 循环参数
! Where结构中也不可调用子程序
! real :: temp
! arr = 0.
! where(arr == 0.)
! call random_seed()
! call random_number(temp)
! arr = temp
! end where
do i = 1, 100*arrayLength
do j = 1, arrayLength
do k = 1, 3*arrayLength
call random_seed()
call random_number(arr(i, j, k)) ! 对数组进行循环赋值
end do
end do
end do
arr = arr*1300 ! 调整数组范围
arrayTemp = arr ! 拷贝数组
write(*, *) 'Do循环判断开始, 数组总和为', sum(arr) ! 参照
do i = 1, 100*arrayLength
do j = 1, arrayLength
do k = 1, 3*arrayLength
if (arr(i, j, k) > 1000.) arr(i, j, k) = 1000. ! DO判断
end do
end do
end do
write(*, *) 'Do循环判断结束, 数组总和为', sum(arr) ! 输出结果
write(*, *) 'Where结构开始' ! 参照
where(arrayTemp > 1000.)
arrayTemp = 1000. ! WHERE判断, Where结构真的方便...
end where
write(*, *) 'Where结构结束, 数组总和为', sum(arrayTemp) ! 输出结果
end subroutine Exercises8_12
! 习题8-13 计算年平均温度
subroutine Exercises8_13
implicit none
real, dimension(6, 6) :: array ! 定义数组
integer :: i ! 循环参数
integer :: errorLevel ! 错误码
character(len=20) :: inforError ! 错误信息
open(unit=1, status='scratch', iostat=errorLevel, iomsg=inforError) ! 打开临时文件
if (errorLevel /= 0) stop inforError ! 判断打开成功与否
write(1, 1) 68.2, 72.1, 72.5, 74.1, 74.4, 74.2, & ! 输入统计数据
&69.4, 71.1, 71.1, 71.9, 73.1, 73.6, &
&68.9, 70.5, 70.9, 71.5, 72.8, 73.0, &
&68.6, 69.9, 70.4, 70.8, 71.5, 72.2, &
&68.1, 69.3, 69.8, 70.2, 70.9, 71.2, &
&68.3, 68.8, 69.6, 70.0, 70.5, 70.9
1 format(6F5.1) ! 设定统计数据输入格式
rewind(unit=1) ! 回到文件开头
read(1,*) (array(i, :), i = 1, 6) ! 给数组赋值
2 format(1X, F4.1, A8, ' 的平均气温为 ', F4.1, ' ℃') ! 设定结果输出格式
do i = 1, 6
write(*, 2) 29.5 + real(i)*0.5, '°N lat ', sum(array(i, :))/6. ! 输出纬度平均气温
end do
do i = 1, 6
write(*, 2) 89.5 + real(i)*0.5, '°W long', sum(array(:, i))/6. ! 输出经度平均气温
end do
write(*, 3) '所有统计地点年平均温度为 ', sum(array)/36., ' ℃' ! 输出所有地区平均气温
3 format(1X, A, F4.1, A) ! 设定结果输出格式
close(unit=1) ! 关闭文件, 此句可忽略
end subroutine Exercises8_13
! 习题8-14 矩阵乘法, 注意, F03虽然可以自动分配大小, 但是还是得先给他分配一个大小
subroutine Exercises8_14
implicit none
! real, dimension(:,:), allocatable :: MatrixMultplication
real, dimension(:, :), allocatable :: array1, array2, array3
integer :: i, j, k
integer :: a, b
integer :: errorLevel ! 错误码
character(len=20) :: inforError ! 错误信息
open(unit=1, status='scratch', iostat=errorLevel, iomsg=inforError) ! 打开临时文件
if (errorLevel /= 0) stop inforError ! 判断打开成功与否
open(unit=2, status='scratch', iostat=errorLevel, iomsg=inforError) ! 打开临时文件
if (errorLevel /= 0) stop inforError ! 判断打开成功与否
! 习题8-14使用
! 匹配的, 输入矩阵1
! write (1, '(A)') " 2, 2",&
! &" 3.0, -1.0",&
! &" 1.0, 2.0"
! ! 不匹配的, 输入矩阵1
! ! write (1, '(A)') " 2, 2, 3",&
! ! &" 3.0, -1.0, 1.0",&
! ! &" 1.0, 2.0, 1.0"
! ! 输入矩阵2
! write (2, '(A)') " 2, 2",&
! &" 1.0, 4.0",&
! &" 2.0, -3.0"
! 习题8-15使用
write (1, '(A)') " 2, 4",&
&" 1.0, -5.0, 4.0, 2.0",&
&" -6.0, -4.0, 2.0, 2.0"
write (2, '(A)') " 4, 3",&
&" 1.0, -2.0, -1.0",&
&" 2.0, 3.0, 4.0",&
&" 0.0, -1.0, 2.0",&
&" 0.0, -3.0, 1.0"
rewind(unit=1) ! 回到文件1第一行
rewind(unit=2) ! 回到文件2第一行
read(1, *) a, b ! 读取矩阵1大小
allocate(array1(a, b), stat=errorLevel) ! 分配矩阵1
if (errorLevel /= 0) print *, "array1: Allocation request denied"
read(1, *) (array1(i, :), i = 1, a) ! 给矩阵1赋值
read(2, *) a, b ! 读取矩阵2大小
allocate(array2(a, b), stat=errorLevel) ! 分配矩阵2
if (errorLevel /= 0) print *, "array2: Allocation request denied"
read(2, *) (array2(i, :), i = 1, a) ! 给矩阵2赋值
! 关闭打开的临时文件(实际上一般可以不管)
close(unit=1)
close(unit=2)
if (size(array1, 2) /= size(array2, 1)) stop "警告! 数组维度不匹配!"! 检测是否可进行矩阵乘法
allocate(array3(size(array1, 1), size(array2, 2)), stat=errorLevel) ! 分配矩阵3
if (errorLevel /= 0) print *, "array3: Allocation request denied"
! array3 = MatrixMultplication(array1, array2)
array3 = 0. ! 初始化矩阵3
! 习题8-14, 8-15使用
! FORALL结构给矩阵3赋值, 注意这里会警告重复计算, 别管他, 设计就是这样
! forall(i = 1:size(array1, 1), j = 1:size(array1, 2), k = 1:size(array2, 2))
! array3(i, k) = array1(i, j)*array2(j, k) + array3(i, k)
! end forall
! 习题8-16 使用(改写后的)
array3 = matmul(array1, array2)
! 关闭数组, 释放资源
if (allocated(array2)) deallocate(array2, stat=errorLevel)
if (errorLevel /= 0) print *, "array2: Deallocation request denied"
if (allocated(array1)) deallocate(array1, stat=errorLevel)
if (errorLevel /= 0) print *, "array1: Deallocation request denied"
write(*, 2) (array3(i, :), i = 1, size(array3, 1)) ! 输出矩阵乘法计算结果
2 format(2(1X, F7.2, 1X))
if (allocated(array3)) deallocate(array3, stat=errorLevel)
if (errorLevel /= 0) print *, "array3: Deallocation request denied"
! contains
! function MatrixMultplication(arrayInput1, arrayInput2)
! implicit none
! real, dimension(:,:), allocatable, intent(in) :: arrayInput1, arrayInput2
! real, dimension(:,:), allocatable :: MatrixMultplication
! integer :: i, j, k
! if (size(arrayInput1, 2) /= size(arrayInput2, 1)) stop "警告! 数组维度不匹配!"
! MatrixMultplication = 0.
! forall(i = 1:size(arrayInput1, 1), j = 1:size(arrayInput1, 2), k = 1:size(arrayInput2, 2))
! MatrixMultplication(i, k) = arrayInput1(i, j)*arrayInput2(j, k) + MatrixMultplication(i, k)
! end forall
! end function MatrixMultplication
end subroutine Exercises8_14
! 习题8-17 相对极大(鞍点)
! 根据鞍点的定义, 可以直接遍历程序中每个点,,, 不过这样好像效率不太高的样子
! 1. 鞍点比上一行的值都大, 比下一行的值都大, 比左一列的值都大, 比右一列的值都大
! 2. 检测筛选之后的值
! 所以, 需要有一个输入矩阵, 有一个标定矩阵, 确定其位置的, i, j, 以及坐标数组(可扩展)
subroutine Exercises8_17
implicit none
real, dimension(:, :), allocatable :: arrayInput ! 输入矩阵
logical, dimension(:, :), allocatable :: arrayCheck ! 检测矩阵
integer :: i, j ! 循环参数
integer, dimension(:, :), allocatable :: arrayPoint ! 鞍点坐标
integer :: sizeX, sizeY ! 矩阵大小
integer :: errorLevel ! 错误码
integer :: num ! 鞍点个数
character(len=20) :: inforError ! 错误信息
! 变量初始化
num = 0
open(unit=1, status='scratch', iostat=errorLevel, iomsg=inforError) ! 打开临时文件
if (errorLevel /= 0) stop inforError ! 确认临时文件成功打开
write (1, '(A)') " 8, 8",& ! 向临时文件写入矩阵
&" 2.0, -1.0, -2.0, 1.0, 3.0, -5.0, 2.0, 1.0",&
&" -2.0, 0.0, -2.5, 5.0, -2.0, 2.0, 1.0, 0.0",&
&" -3.0, -3.0, -3.0, 3.0, 0.0, 0.0, -1.0, -2.0",&
&" -4.5, -4.0, -7.0, 6.0, 1.0, -3.0, 0.0, 5.0",&
&" -3.5, -3.0, -5.0, 0.0, 4.0, 17.0, 11.0, 5.0",&
&" -9.0, -6.0, -5.0, -3.0, 1.0, 2.0, 0.0, 0.5",&
&" -7.0, -4.0, -5.0, -3.0, 2.0, 4.0, 3.0, -1.0",&
&" -6.0, -5.0, -5.0, -2.0, 0.0, 1.0, 2.0, 5.0"
rewind(unit=1) ! 返回初始行
read(1, *) sizeY, sizeX ! 获得矩阵大小
allocate(arrayInput(sizeX, sizeY), stat=errorLevel) ! 定义输入矩阵
if (errorLevel /= 0) print *, "arrayInput: Allocation request denied"
allocate(arrayCheck(sizeX, sizeY), stat=errorLevel) ! 定义检测矩阵
if (errorLevel /= 0) print *, "arrayCheck: Allocation request denied"
allocate(arrayPoint(2, 1), stat=errorLevel) ! 定义鞍点坐标数组
if (errorLevel /= 0) print *, "arrayPoint: Allocation request denied"
read(1, *) (arrayInput(i, :), i = 1, sizeY) ! 给输入矩阵赋值
close(unit=1, iostat=errorLevel, iomsg=inforError) ! 关闭临时文件
if (errorLevel /= 0) stop inforError ! 检测是否关闭成功
arrayCheck = .FALSE. ! 初始化检测矩阵
do i = 2, sizeX - 1 ! 标记行间最大值
arrayCheck(i, 2:sizeY-1) = (arrayInput(i, 2:sizeY-1) > arrayInput(i-1, 2:sizeY-1)) .and. &
&(arrayInput(i, 2:sizeY-1) > arrayInput(i+1, 2:sizeY-1))
end do
do j = 2, sizeY - 1 ! 标记行列间最大值
arrayCheck(2:sizeX-1, j) = (arrayInput(2:sizeX-1, j) > arrayInput(2:sizeX-1, j-1)) .and. &
&(arrayInput(2:sizeX-1, j) > arrayInput(2:sizeX-1, j+1)) .and. &
& arrayCheck(2:sizeX-1, j)
end do
forall (i=2:sizeX-1, j=2:sizeY-1, arrayCheck(i, j)) ! 标记点与左上比较并迭代
arrayCheck(i, j) = (arrayInput(i, j) > arrayInput(i-1, j-1))
end forall
forall (i=2:sizeX-1, j=2:sizeY-1, arrayCheck(i, j)) ! 标记点与左下比较并迭代
arrayCheck(i, j) = (arrayInput(i, j) > arrayInput(i+1, j-1))
end forall
forall (i=2:sizeX-1, j=2:sizeY-1, arrayCheck(i, j)) ! 标记点与右下比较并迭代
arrayCheck(i, j) = (arrayInput(i, j) > arrayInput(i+1, j+1))
end forall
forall (i=2:sizeX-1, j=2:sizeY-1, arrayCheck(i, j)) ! 标记点与右上得到鞍点
arrayCheck(i, j) = (arrayInput(i, j) > arrayInput(i-1, j+1))
end forall
if (allocated(arrayInput)) deallocate(arrayInput, stat=errorLevel) ! 关闭输入矩阵
if (errorLevel /= 0) print *, "arrayInput: Deallocation request denied"
do i = 2, sizeX - 1
do j = 2, sizeY - 1
if (arrayCheck(j, i) .eqv. .TRUE.) then
num = num + 1
arrayPoint(:, num) = (/i, j/) ! 遍历检测矩阵, 得到鞍点坐标
end if
end do
end do
if (allocated(arrayCheck)) deallocate(arrayCheck, stat=errorLevel) ! 关闭检测矩阵
if (errorLevel /= 0) print *, "arrayCheck: Deallocation request denied"
write(*, 1) '矩阵中的鞍点坐标如下:', (arrayPoint(:, i), i = 1, num) ! 输出结果
1 format(1X, A/&
&6('(', I2, ',', I2, ') ')) ! 定义输出格式
if (allocated(arrayPoint)) deallocate(arrayPoint, stat=errorLevel) ! 关闭坐标数组
if (errorLevel /= 0) print *, "arrayPoint: Deallocation request denied"
end subroutine Exercises8_17
! 习题8-18 金属盘温度分布, 总感觉这题好像有点难度(如果要优化算法的话, 直接穷举确实容易倒是, 就是浪费计算资源, 模仿自然温度扩散写个循环吧)
! 基本路线原则如下:
! 1. 计算前进方向一格周围的格子情况, 以确定下一次前进方向, 并累计前进次数
! 1.1 当前进次数为0时, 朝方向1前进1格, 前进次数+1,
! 1.2 当前进次数大于1时:
! 1.2.1 检测前进方向格子状态
! 1.2.1.1 如果为0(未计算), 前进一格, 并检测 优先侧方向 格子状态
! 1.2.1.1.1 如果为0(未计算), 则下次前进方向为 优先侧方向
! 1.2.1.1.2 如果为1(已计算)或-1(固定值), 检测 前进方向 格子状态
! 1.2.1.1.2.1 如果为0(未计算), 则下次前进方向为 当前方向
! 1.2.1.1.2.2 如果为1(已计算)或-1(固定值). 检测 非优先侧方向
! 1.2.1.1.2.2.1 如果为0(未计算), 则下次前进方向为 非优先侧方向
! 1.2.1.1.2.2.2 如果为1(已计算)或-1(固定值), 这时候无路可走, 检测所有格子状态
! 1.2.1.1.2.2.2.1 所有格子为1, 检测此次迭代计算结果
! 1.2.1.1.2.2.2.1.1 符合条件, 退出
! 1.2.1.1.2.2.2.1.1 不符合条件, 初始化检测状态, 进行下一次循环
! 1.2.1.1.2.2.2.2 存在格子不为1, 选定可选位置中
! 1.2.1.2 如果为1(已计算或固定值)
! 突然发现越来越复杂了, 舍弃这种思路
! 按照距离优先计算, 优先计算离热源近的方块, 但问题是这样只能适应单热源的部分, 无法计算多热源(说得好像之前那种就能计算多热源了)
! 为了处理多热源的问题, 引入时序步进的想法, 假定热量传递速度相同
! 这样就是多热源画同样的同心圆, 在此范围内未被计算过的方块会被标记为2
! 对标记为2的方块进行统计, 周围被标记为0的格子数为2的, 被标记为3
! 优先计算标记为2的方块, 再计算标记为3的方块
! 步进计算, 直到无方块被标记
! 从距离的角度, 选取欧式距离似乎最为合适, 选取曼哈顿距离则计算最为简单, 最终还是选取了欧式距离
! 有时候感觉自己就是个ZZ, 直接算出来每个点的距离, 然后排序计算不就好了, 干嘛费那么大劲还去考虑步进每次包括的点位和点位计算顺序
subroutine Exercises8_18
implicit none
real, dimension(10, 10) :: metalPlate, metalPlateOld ! 金属盘表面温度, 及对照组
integer, dimension(63, 2) :: pointList ! 优化后的计算顺序
integer :: i, j, k ! 循环参数(什么都干的那种)
real :: temp, distense ! 用于判断是否符合要求
! 初始化计算顺序用数组
metalPlate = 0.
metalPlate(2:9, 2:9) = 1.
metalPlate(3, 8) = 0.
! 计算需要计算温度处距离热源的距离
forall (i = 2: 9, j = 2: 9, metalPlate(i, j) == 1.)
metalPlate(i, j) = (i-3)**2 + (j-8)**2
end forall
i = 0
! 输出计算顺序
do
i = i + 1
pointList(i, :) = minloc(metalPlate, mask=metalPlate>0.)
metalPlate(pointList(i, 1), pointList(i, 2)) = 0.
if (all(metalPlate == 0.)) exit
end do
! 初始化金属盘表面温度
metalPlate = 20.
metalPlate(2:9, 2:9) = 50.
metalPlate(3, 8) = 100.
j = 0
! 开始迭代计算金属盘表面温度
do
distense = 0.
do i = 1, size(pointList, 1)
temp = metalPlate(pointList(i, 1), pointList(i, 2))
metalPlate(pointList(i, 1), pointList(i, 2)) = (metalPlate(pointList(i, 1)+1, pointList(i, 2)) + &
&metalPlate(pointList(i, 1), pointList(i, 2)+1) + metalPlate(pointList(i, 1)-1, pointList(i, 2)) &
&+ metalPlate(pointList(i, 1), pointList(i, 2)-1))/4.
if (abs(temp-metalPlate(pointList(i, 1), pointList(i, 2))) > distense)&
& distense = abs(temp-metalPlate(pointList(i, 1), pointList(i, 2)))
end do
if (distense < 0.01) exit
j = j + 1
end do
write(*, '(" 共经历", I3, "次迭代计算")') j
write(*, *) '节点(5, 5)处稳定状态下温度为', metalPlate(5, 5)
write(*, '(" ", 10F7.2)') (metalPlate(i, :), i = 1, 10)
! 初始化金属盘表明温度
metalPlate = 20.
metalPlate(2:9, 2:9) = 50.
metalPlate(3, 8) = 100.
k = 0
! 开始迭代计算金属盘表面温度
do
metalPlateOld = metalPlate
forall(i = 1:10, j = 1:10, metalPlate(i, j) /= 20. .and. metalPlate(i, j) /= 100.)
metalPlate(i, j) = (metalPlate(i+1, j) + metalPlate(i, j+1) + metalPlate(i-1, j) + metalPlate(i, j-1))/4.
end forall
k = k + 1
if (all(abs(metalPlate-metalPlateOld) < 0.01)) exit
end do
write(*, '(" 共经历", I3, "次迭代计算")') k
write(*, *) '节点(5, 5)处稳定状态下温度为', metalPlate(5, 5)
write(*, '(" ", 10F7.2)') (metalPlate(i, :), i = 1, 10)
end subroutine Exercises8_18