Fortran入门.07.判断,副程序 - 北方连萌

Fortran入门.07.判断,副程序

If

写法

if (条件) then
  满足条件时执行的部分
else
  不满足条件时执行的部分
end if

不满足条件时不执行,则

if (条件) then
  满足条件时执行的部分
end if

再精简

if (条件) 满足条件时执行的部分

判断条件增加

if (条件1) then
  执行部分1
else if (条件2) then
  执行部分2
else if (条件3) then
  执行部分3
else
  以上条件皆不满足时执行的部分
end if

条件式中不等式的不等号写法

写法意味
<小于
<=小于等于
==等于
>=大于等于
>大于
/=不等于

用例

素数筛(埃拉托斯特尼筛法)

找出[2, n]之间的素数
第一步:写n个0
第二步:将2的倍数(除2自身)全部标记为1
第三步:将3的倍数(除3自身)全部标记为1
......
第n步:将n的倍数(除n自身)全部标记为1
剩下还是0的即为素数

program abc
  implicit none
  integer :: i, j, k = 1
  integer, parameter :: n = 2000
  integer, dimension(n) :: a, b

  a = 0
  do i = 2, n
    if (a(i) == 0) then
      b(k) = i
      k = k + 1
      do j = i * i, n, i
        a(j) = 1
      end do
    end if
  end do
  print '(1x, 15i6)', b(1: k - 1)

  print *,'press any key to continue'
  read *,
    
  stop
    
end
      2     3     5     7    11    13    17    19    23    29    31    37    41    43    47
     53    59    61    67    71    73    79    83    89    97   101   103   107   109   113
    127   131   137   139   149   151   157   163   167   173   179   181   191   193   197
    199   211   223   227   229   233   239   241   251   257   263   269   271   277   281
    283   293   307   311   313   317   331   337   347   349   353   359   367   373   379
    383   389   397   401   409   419   421   431   433   439   443   449   457   461   463
    467   479   487   491   499   503   509   521   523   541   547   557   563   569   571
    577   587   593   599   601   607   613   617   619   631   641   643   647   653   659
    661   673   677   683   691   701   709   719   727   733   739   743   751   757   761
    769   773   787   797   809   811   821   823   827   829   839   853   857   859   863
    877   881   883   887   907   911   919   929   937   941   947   953   967   971   977
    983   991   997  1009  1013  1019  1021  1031  1033  1039  1049  1051  1061  1063  1069
   1087  1091  1093  1097  1103  1109  1117  1123  1129  1151  1153  1163  1171  1181  1187
   1193  1201  1213  1217  1223  1229  1231  1237  1249  1259  1277  1279  1283  1289  1291
   1297  1301  1303  1307  1319  1321  1327  1361  1367  1373  1381  1399  1409  1423  1427
   1429  1433  1439  1447  1451  1453  1459  1471  1481  1483  1487  1489  1493  1499  1511
   1523  1531  1543  1549  1553  1559  1567  1571  1579  1583  1597  1601  1607  1609  1613
   1619  1621  1627  1637  1657  1663  1667  1669  1693  1697  1699  1709  1721  1723  1733
   1741  1747  1753  1759  1777  1783  1787  1789  1801  1811  1823  1831  1847  1861  1867
   1871  1873  1877  1879  1889  1901  1907  1913  1931  1933  1949  1951  1973  1979  1987
   1993  1997  1999
 press any key to continue

利用迭代法求数组的最大 / 小值

先将数组中的第一个数同时指定为最大值和最小值
对第2个数,如果其大 / 小于前次指定的最大 / 小值,则将其指定为新最大 / 小值
......
对第n个数,如果其大 / 小于前次指定的最大 / 小值,则将其指定为新最大 / 小值
比较完成,输出最新一轮比较出的最大 / 小值

program abc
  implicit none
  integer :: i, n
  integer, parameter :: nn = 2000
  real, dimension(nn) :: a
  real :: maax, miin

  print *, 'n = '
  read *, n
  do i = 1, n
    print '(1x, 2ha(, i1, 2h)= )', i
    read *, a(i)
  end do

  maax = a(1)
  miin = a(1)
  do i = 2, n
    if(a(i) > maax) maax = a(i)
    if(a(i) < miin) miin = a(i)
  end do

  print *, 'maximum', maax
  print *, 'minimum', miin

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 n =
6
 a(1)=
54515.231
 a(2)=
45345.11
 a(3)=
5143515012
 a(4)=
1231.011
 a(5)=
1513.01515
 a(6)=
1631194.001144
 maximum   5.14351514E+09
 minimum   1231.01099
 press any key to continue

交换法(冒泡排序)

对第1到第n - 1项
将每一项与后项比较,如果该项较后项大,则将两项交换位置
此时最后一项应是最大的,故对第1到第n - 2 项
将每一项与后项比较,如果该项较后项大,则将两项交换位置
......
对第1项
与后项比较,大则换位
排序完成

program abc
  implicit none
  integer :: i, n, ter
  integer, parameter :: nn = 2000
  real, dimension(nn) :: a
  real :: bigger

  print *, 'n = '
  read *, n
  print *, 'a = '
  read *, a(1: n)

  do ter = n - 1, 1, -1
    do i = 1, ter
      if (a(i) > a(i + 1)) then
        bigger= a(i)
        a(i) = a(i + 1)
        a(i + 1) = bigger
      end if
    end do
  end do

  print *, a(1: n)

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 n =
6
 a =
5115
6161
52152
121
3513
1315
   121.000000       1315.00000       3513.00000       5115.00000       6161.00000       52152.0000
 press any key to continue

JUMP

Exit

出现于循环节(do)之中,用来退出循环

do 循环条件
  if (循环执行到某次总会满足的一个条件) then
    循环执行的部分
    exit
  end if
end do

多重循环的场合

命名1: do
  命名2: do
          if (...) exit 命名1          
          if (...) exit 命名2
  end do
end do

不加指定退出哪一个循环节则只退出最内侧的那个

Go to

跳到指定序号的那一行

   if (...) go to 10
     ...
10 print *, '10'

Cycle

出现于循环节(do)之中,用来初始化循环

do ...
  if (...) cycle
  循环执行的部分
end do

Do while

写法

满足条件时一直重复

do while (条件)
  循环执行的部分
end do

用例

二分查找算法(整数)

确定搜索的上下限
求上下限的平均数并取整
如果要查找的数正好是平均数,则过程结束
如果要查找的数在平均数与上 / 下限之间,则将平均数指定为新的最下 / 上限
......
由于上限 - 下限 = 2时再操作一次平均数必定是最终结果
故循环条件是上限 > 下限 + 1

例如求分段函数的取值

$$\begin{array}\ y=\begin{cases}300\quad x\in (0,2]\\\ 550\quad\ x\in \ (2,3]\\\ 880\quad\ x\in \left(3,5\right]\\\ 1100\quad\ x\in \left(5,7\right]\\\ 1600+500\times \left(\frac{x-7.000001}{5}\right)\quad\ x\in \left(7,\infty \right)\ \end{cases}\end{array}\begin{array}\ x\in N*\end{array}$$


program abc
  implicit none
  integer :: d, u, m, t
  integer, parameter :: n = 4
  integer, parameter :: q = 3
  integer, dimension(n) :: x, y
  integer, dimension(q) :: a, b

  x = (/2, 3, 5, 7/)
  y = (/300, 550, 880, 1100/)
  
  do t = 1, 3
    print '(1x, a2, i2, a2)', 'x(', t, ')='
    read *, a(t)

    if(a(t) <= x(1)) then
      b(t) = y(1)
    else if(a(t) > x(n)) then
      b(t) = 1600 + 500 * int((a(t) - 7.000001) / 5)
    else
      d = 1
      u = n
      do while(d + 1 < u)
        m = (u + d) / 2
        if (a(t) == x(m)) then
          u = m
          exit
        else if(a(t) < x(m)) then
          u = m
        else
          d = m
        end if
      end do
      b(t) = y(u)
    end if
  end do

  do t = 1, 3
    print '(1x, a2, i2, a2, i5)', 'y(', t, ')=', b(t)
  end do

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 x( 1)=
1
 x( 2)=
6
 x( 3)=
15
 y( 1)=  300
 y( 2)= 1100
 y( 3)= 2100
 press any key to continue

当然,将输入的值与有界部分每段的上界比较,然后找出该值所在的区间的方法也行得通

program abc
  implicit none
  integer :: i, t
  integer, parameter :: n = 4
  integer, parameter :: q = 3
  integer, dimension(n) :: x, y
  integer, dimension(q) :: a, b

  x = (/2, 3, 5, 7/)
  y = (/300, 550, 880, 1100/)
  
  do t = 1, 3
    print '(1x, a2, i2, a2)', 'x(', t, ')='
    read *, a(t)

    if(a(t) > x(n)) then
      b(t) = 1600 + 500 * int((a(t) - 7.000001) / 5)
    else
      do i = 1, 4
        if(a(t) <= x(i)) then
          b(t) = y(i)
          exit
        end if
      end do
    end if
  end do

  do t = 1, 3
    print '(1x, a2, i2, a2, i5)', 'y(', t, ')=', b(t)
  end do

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 x( 1)=
1
 x( 2)=
6
 x( 3)=
13
 y( 1)=  300
 y( 2)= 1100
 y( 3)= 2100
 press any key to continue

但是设想极端情况,比如要查找的数在有界部分的很后端,那么用逐个对比的方法就会比不断二分来得慢了

副程序

程序较长时,分主次更合适。
副程序一般用来摞subroutine(子程序)和function(函数)。
副程序是独立于主程序的,所以咱一般把它写在

end program

的后一行。
也是因此,需要在主程序的宣言中宣言副程序要用到的各参数,以便调动副程序中的参数到主程序中。

如果是想要在副程序中再创建副程序,步骤与上述相同,也就是写在

end subroutine

或者

end function

的后行。

你会注意到这个时候同一个宣言要在主副里头都要写一道,很麻烦。如果想要省略一道,可以用

contains

把external(外部的)subroutine或function变成internal(内部的)subroutine或function。

具体参照以下内容。

写法

Subroutine / function ( external )

在主程序的宣言后追加

interface
  subroutine 子程序名(用到的数)
  implicit none
  型, intent (in) :: 输入用数扔这
  型, intent (out) :: 输出用数扔这
  型, intent (inout) :: 输入输出两用数扔这
  end subroutine 子程序名
end interface
interface
  function 函数名(用到的数)
  implicit none
  型, intent (in) :: 输入用数扔这
  型, intent (out) :: 输出用数扔这
  型, intent (inout) :: 输入输出两用数扔这
  end function 函数名
end interface

在主程序结束以后写

subroutine 子程序名(用到的数)
implicit none
型, intent (in) :: 列出输入用数
型, intent (out) :: 列出输出用数
型, intent (inout) :: 列出输入输出两用数
对subroutine中要用到的其他进行型宣言
  处理
  结果代入输出用数或输入输出两用数
  return
end subroutine 子程序名
function 函数名(用到的数)
implicit none
型, intent (in) :: 列出输入用数
型, intent (out) :: 列出输出用数
型, intent (inout) :: 列出输入输出两用数
对function中要用到的其他进行型宣言
  处理
  函数名 = 处理结果
  return
end function 函数名

Subroutine / function ( internal )

主程序stop以后用写包含的若干副程序,然后end掉主程序

program 主程序名
  implicit none
  型宣言
  处理
  stop

  contains !包含

  subroutine 子程序名(用到的数)
    implicit none
    型, intent (in) :: 列出输入用数
    型, intent (out) :: 列出输出用数
    型, intent (inout) :: 列出输入输出两用数
    对subroutine中要用到的其他进行型宣言
    处理
    结果代入输出用数或输入输出两用数
    return
  end subroutine 子程序名

  function 函数名(用到的数)
  implicit none
  型, intent (in) :: 列出输入用数
  型, intent (out) :: 列出输出用数
  型, intent (inout) :: 列出输入输出两用数
  对function中要用到的其他进行型宣言
    处理
    函数名 = 处理结果
    return
  end function 函数名

end program

用例

分数的加法

呼出子程序时,用以下表达

call 子程序名(用到的数)

如计算

$$\frac{a}{b} + \frac{c}{d} = \frac{e}{f}$$

program abc
  implicit none
  integer :: a, b, c, d, e, f
  interface
    subroutine fs(a, b, c, d, e, f) !外部子程序与主程序相互独立,可以重复使用主程序中已经使用到的字母
    implicit none
    integer, intent(in) :: a, b, c, d
    integer, intent(out) :: e, f
    end subroutine
  end interface

  print *, 'a, b, c, d = '
  read *, a, b, c, d
  call fs(a, b, c, d, e, f) !呼出子程序
  print *, 'e = ', e
  print *, 'f = ', f

  print *,'press any key to continue'
  read *,
    
  stop
    
end
!以下是通分运算再约分的子程序
subroutine fs(a, b, c, d, e, f)
  implicit none
  integer, intent(in) :: a, b, c, d
  integer, intent(out) :: e, f
  integer :: g
  interface
    subroutine gcd(a, b, c)
      implicit none
      integer, intent(in) :: a, b
      integer, intent(out) :: c
    end subroutine
  end interface

  e = a * d + b * c
  f = b * d
  call gcd(e, f, g) !子程序中呼出子程序
  e = e / g
  f = f / g
  return

end subroutine
!以下是辗转相除法求最大公约数的子程序
subroutine gcd(a, b, c)
  implicit none
  integer, intent(in) :: a, b
  integer, intent(out) :: c
  integer :: d, e, f

  d = a
  e = b
  do
    f = mod(d, e)
    if(f == 0) then
      c = e
      exit
    else
      d = e
      e = f
    end if
  end do
  return

end subroutine
 a, b, c, d =
64641,65465,622,8485
 e =    117839623
 f =    111094105
 press any key to continue
向量的夹角

根据函数的定义,函数名在此是唯一的输出值,故在宣言时只用intent(in),函数名不用宣言为intent(out)
另外函数的呼出也相对简单,形如笔算的f(x)型就可以

如计算3维向量a, b的夹角

program abc
  implicit none
  real, dimension(3) :: a, b
  real :: s, t
  interface
    function dp(a, b, n)
      implicit none
      real :: dp
      real, dimension(*), intent(in) :: a, b
      integer, intent(in) :: n
    end function
  end interface

  print *, 'Enter a 3 dimensional vector a'
  read *, a
  print *, 'Enter a 3 dimensional vector b'
  read *, b
  s = dp(a, b, 3) / sqrt(dp(a, a, 3) * dp(b, b, 3)) !直接把定义好的函数代入运算
  t = acos(s) * 180 / 3.1415926535
  print *, 'The angle between a and b is', t, 'degree(s).'

  print *,'press any key to continue'
  read *,
    
  stop
    
end
!以下时n维向量的内积计算函数
function dp(a, b, n)
  implicit none
  real :: dp
  real, dimension(*), intent(in) :: a, b
  integer, intent(in) :: n
  integer :: i
  real :: s

  s = 0.0
  do i = 1, n
    s = a(i) * b(i) + s
  end do
  dp = s
  return

end function
 Enter a 3 dimensional vector a
44,55,66
 Enter a 3 dimensional vector b
11,77,99
 The angle between a and b is   22.1832886     degree(s).
 press any key to continue
二项式系数
二分法求根

确定初始的x1,x2,使f(x1)<0,f(x2)>0
求x1,x2的平均数x
若f(x)<0,则将x作为新的x1,否则将其作为新的x2重复操作
......
达到指定回数时,停止操作

如求

$$\cos\ x-x=0$$

的根

program abc
  implicit none
  real :: x, y, x1, x2
  integer :: r

  x1 = 0.0
  x2 = 1.0
  do r = 1, 10000
    x = (x1 + x2) / 2.0
    y = f(x)
    if(y < 0) then
      x1 = x
    else
      x2 = x
    end if
  end do

  print *, 'x = ', x

  print *,'press any key to continue'
  read *,
    
  stop
!主程序中包含以下函数
  contains
  function f(x)
    implicit none
    real :: f
    real, intent(in) :: x

    f = x - cos(x)
    return

  end function
    
end
 x =   0.739085078
 press any key to continue
牛顿法求根

先看一个简单的例子

$$x^{2}-99=0$$

这就相当于是求面积是99的正方形的边长

于是可以有以下思路

首先找到一组长宽,使其乘积等于给定的正方形面积
接下来把长宽的平均数指定为新长,面积除以新长指定为新宽
......
直到长宽之差几乎为0,结束过程

program abc
  implicit none
  real :: l, w, x2
  real, parameter :: p = 0.0000001

  print *, 'square = '
  read *, x2

  l = 1.0
  w = x2
  do while(abs((l - w) / l) > p)
    l = (l + w) * 0.5
    w = x2 / l
  end do

  print *, 'root = ', l

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 square =
99
 root =    9.94987488
 press any key to continue

如果要求任意函数与x轴的交点,如

$$\cos\ x-x=0$$

则使用牛顿法

通过不断作切线,直到f(x)几乎为0

$$x_{n+1} = x_{n} - \frac{f(x_{n})}{f^{\prime}(x_{n})}$$

逼近思想与与上述长方化正方是一致的

数值求导可以利用导数的定义求,为避免麻烦此处直接使用笔算的表达式

program abc
  implicit none
  real :: x, y, dy, p = 0.000001
  integer :: r

  x = 1.0
  do r = 1, 100
    y = f(x)
    dy = df(x)
    if(abs(dy) < p) then
      print *, 'Stationary Point'
      stop
    end if
    x = x - y / dy
    if(abs(y) < p) exit !y值足够小
    if(abs(y / dy) < p) exit !x(n+1)和x(n)的距离足够小
  end do

  print *, 'x = ', x
  
  print *,'press any key to continue'
  read *,
    
  stop

  contains

  function f(x)
    implicit none
    real :: f
    real, intent(in) :: x

    f = cos(x) - x
    return

  end function

  function df(x)
    implicit none
    real :: df
    real, intent(in) :: x

    df = -sin(x) - 1
    return

  end function
    
end
 x =   0.739085138
 press any key to continue
辛普森积分法

$$\int_{a}^{b} f(x) d x \approx \frac{h}{3}\left[f\left(x_{0}\right)+2 \sum_{j=1}^{n / 2-1} f\left(x_{2 j}\right)+4 \sum_{j=1}^{n / 2} f\left(x_{2 j-1}\right)+f\left(x_{n}\right)\right]$$
其中
$$h=(b-a) / n \\x_{i}=a+i h(i=0, \ldots, n)$$

program abc
  implicit none
  integer :: j, n
  real :: a, b, h, s1, s2, s
  interface
    function f(x)
      implicit none
      real :: f
      real, intent(in) :: x
    end function
  end interface

  print *, 'a, b = '
  read *, a, b
  print *, 'n = '
  read *, n

  h = (b - a) / real(n)

  s1 = 0.0
  do j = 2, n - 2, 2
    s1 = s1 + f(a + j * h)
  end do

  s2 = 0.0
  do j = 1, n - 1, 2
    s2 = s2 + f(a + j * h)
  end do

  s = (f(a) + 2.0 * s1 + 4.0 * s2 + f(b)) * h / 3.0

  print *, 's = ', s
  
  print *,'press any key to continue'
  read *,
    
  stop
  
end

function f(x)
  implicit none
  real :: f
  real, intent(in) :: x
  
  f = sin(x)**cos(x)
  return

end function
 a, b =
0,1
 n =
1000000
 s =   0.501324713
 press any key to continue

添加新评论

电子邮件地址不会被公开,评论内容可能需要管理员审核后显示。