if 语句确定稳态

Question

我下面的代码正确求解了函数 u(x,t) 的一维热方程。我现在想找到稳态解，不再随时间变化的解，所以它应该满足u(t+1)-u(t) = 0。找到稳态解的最有效方法是什么？我在下面展示了三种不同的尝试，但我不确定两者是否真的在做我想做的事情。第一个和第三个具有正确的语法，第二个方法由于 if 语句而存在语法错误。由于if结构的变化，每个方法都不一样。

方法一：

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: this checks the solutions at every space point  which I don't think is correct.
    do i = 1,n-1
       if ( v(i) - u(i) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if
    end do

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1

方法二：

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: (This gives an error message since the if statement doesn't have a logical scalar expression, but I want to compare the full arrays v and u as shown.
       if ( v - u .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1

方法三：

program parabolic1
  integer, parameter :: n = 10, m = 20 
  real, parameter :: h = 0.1, k = 0.005 !step sizes
  real, dimension (0:n) :: u,v
  integer:: i,j
  real::pi,pi2
  
  u(0) = 0.0; v(0) = 0.0; u(n) = 0.0; v(n) =0.0
  pi = 4.0*atan(1.0)    
  pi2 = pi*pi 
  
  do i=1, n-1
     u(i) = sin( pi*real(i)*h)
  end do
  
  do j = 1,m 

    do i = 1, n-1
       v(i) = 0.5*(u(i-1)+u(i+1))
    end do 
    t = real(j)*k  !increment in time, now check for steady-state

    !steady-state check: Perhaps this is the correct expression I want to use
       if( norm2(v) - norm2(u) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if

    do i = 1, n-1 !updating solution
       u(i) = v(i)       
    end do 

  end do 
end program parabolic1

Answer 1

无需讨论确定“接近度”的最佳或正确方法（不是真正的编程问题），我们可以专注于方法的 Fortran 部分在做什么。

方法 1 和方法 2 是类似的想法（但在执行时有问题），而方法 3 不同（并且在其他方面有问题）。

还请注意，一般来说，人们想要比较差异的大小 abs(v-u) 而不是（有符号的）差异 v-u。随着迭代的非单调变化，这些是完全不同的。

方法三使用norm2(v) - norm2(u)来检验数组u和v是否相似。这是不正确的。考虑

norm2([1.,0.])-norm2([0.,1.])

而不是更正确的

norm2([1.,0.]-[0.,1.])

方法 2

if ( v - u .LT. 1.0e-7 ) then

存在无效数组表达式的问题，但是“所有点都接近吗？”可以适当地写成

if ( ALL( v - u .LT. 1.0e-7 )) then

（您会在此处找到有关此类数组缩减的其他问题）。

方法 1 尝试了类似的方法，但不正确：

    do i = 1,n-1
       if ( v(i) - u(i) .LT. 1.0e-7 ) then
       print*, 'steady-state condition reached'
       exit 
       end if
    end do

这是不正确的，一方面是错误的，另一方面是微妙的。

首先，当条件第一次测试为真时退出循环，并显示一条消息说已达到稳定状态。这是不正确的：您需要关闭 all 值，而这是测试 any 值关闭。

其次，当条件满足时，你exit。但是你没有 exit 时间迭代循环，你退出了接近测试循环。（exit 没有构造名称离开最里面的 do 构造）。您将处于完全相同的情况，运行在这个最内层的构造之后立即再次满足测试条件是否曾经满足或从未满足（如果曾经满足，您也会收到消息）。您将需要在时间循环中使用构造名称。

我不会展示如何做到这一点（同样这里还有其他问题），因为您还需要修复测试条件，到那时您最好使用 if(all(... （更正方法 2）没有额外的 do 构造。

对于方法 1 和 2，您将得到类似的东西：

if (all(v-u .lt 1e-7)) then
  print *, "Converged"
  exit
end if

对于方法 3：

if (norm2(v-u) .lt. 1e-7) then
  print *, "Converged"
  exit
end if

if 语句确定稳态

if statement to determine steady-state

fortran