[Fortran] 纯文本查看 复制代码
module shared_vars
implicit none
integer, parameter :: dp = kind(0.0d0)
! 全局变量
integer, parameter :: step = 200 ! 每周期步数
integer, parameter :: transient_period = 500 ! 瞬态仿真周期 (稳定后提前终止)
integer, parameter :: stable_period = 10 ! 稳态仿真周期
! 预先分配稳态仿真所需的位移数组
integer, parameter :: steady_state_steps = stable_period * step
real(dp), dimension(steady_state_steps) :: x_values
! 外部激励频率和模型参数
real(dp), parameter :: alpha = 0.4_dp, beta = 0.5_dp, delta = -4.0_dp, F = 3.0_dp
real(dp) :: p
end module shared_vars
program duffingfreqanalysis
include 'link_fnl_shared.h'
use shared_vars
use ivmrk_int
implicit none
integer, parameter :: n = 2 ! 方程数量 (Duffing模型中为2个方程)
real(dp) :: t, y(n)
real(dp) :: period ! 周期
real(dp) :: max_displacement
real(dp) :: p_start, p_end, dp_p ! 频率p的范围和步长
real(dp) :: start_time, end_time, elapsed_time ! 记录时间的变量
external :: fcn ! 外部函数
integer :: i, num_steps
! 输出文件的文件名
open(unit=10, file='amplitude_vs_frequency.txt', status='unknown')
! 设置频率p的范围和步长
p_start = 0.5_dp
p_end = 4.0_dp
dp_p = 0.01_dp
num_steps = int((p_end - p_start) / dp_p) + 1
! 开始计时
call cpu_time(start_time)
!$omp parallel do private(i, p, period, t, y) shared(p_start, dp_p, num_steps)
do i = 1, num_steps
p = p_start + (i - 1) * dp_p
! 计算当前频率对应的周期
period = 2 * 3.141592653589793_dp / p
! 调用瞬态计算子程序
call transient_simulation(t, y, period)
! 调用稳态计算子程序
call steady_state_simulation(t, y, period, max_displacement)
! 输出初始结果到文件
!$omp critical
write (10, '(2F12.6)') p, max_displacement
!$omp end critical
end do
!$omp end parallel do
! 关闭文件
close (unit=10)
! 结束计时
call cpu_time(end_time)
! 计算运行时间
elapsed_time = end_time - start_time
! 输出运行时间到控制台
print *, '运行时间 (秒):', elapsed_time
contains
subroutine transient_simulation(t, y, period)
use shared_vars
implicit none
real(dp), intent(in) :: period
real(dp) :: t, tend, y(n), yprime(n)
integer :: ido, istep
! 初始条件设定
t = 0.0
tend = 0.0
y(1) = 0.0 ! 初始位移
y(2) = 0.0 ! 初始速度
! 初始调用ivmrk
ido = 1
! 瞬态仿真循环
do istep = 1, transient_period
! 更新 tend
tend = tend + period ! 更新结束时间
! 调用ivmrk函数进行积分
call ivmrk(ido, fcn, t, tend, y, yprime)
end do
! 释放空间
ido = 3
call ivmrk(ido, fcn, t, tend, y, yprime)
end subroutine transient_simulation
subroutine steady_state_simulation(t, y, period, max_displacement)
use shared_vars
implicit none
real(dp), intent(in) :: period
real(dp) :: t, tend, y(n), yprime(n), h
integer :: ido, istep
real(dp), intent(out) :: max_displacement
! 初始条件设定
tend = t
h = period / step
! 初始调用ivmrk
ido = 1
! 循环调用ivmrk，计算每个时间步的结果
do istep = 1, steady_state_steps
! 更新 tend
tend = tend + h ! 更新当前时间
! 调用ivmrk函数进行积分
call ivmrk(ido, fcn, t, tend, y, yprime)
! 存储位移到数组
x_values(istep) = y(1) ! y(1)表示位移
end do
! 释放空间
ido = 3
call ivmrk(ido, fcn, t, tend, y, yprime)
! 计算最大位移的绝对值
max_displacement = maxval(abs(x_values))
end subroutine steady_state_simulation
end program duffingfreqanalysis
subroutine fcn(n, t, y, yprime)
use shared_vars
implicit none
integer, intent(in) :: n
real(dp), intent(in) :: t
real(dp), intent(in) :: y(n)
real(dp), intent(out) :: yprime(n)
! 计算 Duffing 方程的导数
yprime(1) = y(2)
yprime(2) = -alpha * y(2) - beta * y(1)**3 + delta * y(1) + F * cos(p * t)
end subroutine fcn