Fortran入门.06.数列,向量,矩阵 - 北方连萌

Fortran入门.06.数列,向量,矩阵

数组

型, dimension(上限: 下限)::命名

斐波那契数列

a0 = 1, a1 = 1, an+1 = an + an-1

program fb
  implicit none
  integer, dimension(0: 10) :: a
  integer :: i

  a(0) = 1.0
  a(1) = 1.0
  do i = 1, 9
    a(i+1) = a(i) + a(i-1)
  end do

  print '(1x,11i4)', (a(i), i = 0, 10)

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

end
    1   1   2   3   5   8  13  21  34  55  89
 press any key to continue

平均值与标准差

program name
  implicit none
  integer :: i, n
  real :: s, x, mu, s2, sigma
  real, dimension(100) :: a

    print *, 'n = '
    read *, n 
    ! 输入n个值
    do i = 1, n
      print '(1x, a2, i2, a2)','a(',i,')=' 
      read *, a(i)
    end do
    ! 计算平均数
    s = 0.0
    do i = 1, n
      s = s + a(i)
    end do
    mu = s / n
    ! 计算标准差
    x=0.0
    do i = 1, n
      x = x + (a(i) - mu)**2
    end do
    s2 = x / (real(n))
    sigma = sqrt(s2)

    print *, 'mean', mu
    print *, 'standard deviation', sigma

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

end program name
 n =
6
 a( 1)=
1513654156
 a( 2)=
514514551
 a( 3)=
513513
 a( 4)=
1235413513
 a( 5)=
61351438979845
 a( 6)=
13212131
 mean   1.02257855E+13
 standard deviation   2.28640885E+13
 press any key to continue

最小二乘法

program abc
  implicit none
  integer :: i, n
  real :: sx, sy, mux, muy, s1, s2, alpha, beta
  integer, dimension(100) :: x, y

  print *, 'n = '
  read *, n 
  ! 对于x输入n个值
  do i = 1, n
    print '(1x, a2, i2, a2)','x(',i,')=' 
    read *, x(i)
  end do
  ! 对于y输入n个值
  do i = 1, n
    print '(1x, a2, i2, a2)','y(',i,')=' 
    read *, y(i)
  end do
  ! 计算x平均数
  sx = 0.0
  do i = 1, n
    sx = sx + x(i)
  end do
  mux = sx / n
  ! 计算y平均数
  sy = 0.0
  do i = 1, n
    sy = sy + y(i)
  end do
  muy = sy / n
  ! 计算系数
  s1 = 0.0
  do i = 1, n
    s1 = s1 + (x(i) - mux) * (y(i) - muy)
  end do
  s2 = 0.0
  do i = 1, n
    s2 = s2 + (x(i) - mux) ** 2
  end do
  beta = s1 / s2
  alpha = muy - beta * mux
  ! 输出结果
  print *, 'mean x = ', mux
  print *, 'mean y = ', muy
  print '(1x, a2, f7.3, a1, f7.3, a1)', 'y=', alpha, '+', beta, 'x'

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 n =
3
 x( 1)=
1
 x( 2)=
2
 x( 3)=
3
 y( 1)=
1
 y( 2)=
3
 y( 3)=
2
 mean x =    2.00000000
 mean y =    2.00000000
 y=  1.000+  0.500x
 press any key to continue

向量

ab同为n维向量时,有

写法意义
a + bai + bi (i = 1, ..., n)
a - bai - bi (i = 1, ..., n)
a * bai * bi (i = 1, ..., n)
a / bai / bi (i = 1, ..., n)
a ^ bai ^ bi (i = 1, ..., n)
dot_product(a, b)点积

内积

program vector
  implicit none
  real, dimension(4) :: a, b, c
  real :: s

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

  c = a * b
  s = dot_product(a, b)

  print '(1x, 5f15.0)', c
  print *, 'dot_product(a, b) =', s

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

end 
 a =
10156,457845,16119,4795
 b =
14545,16197,1616,4641
      147719024.    7415715328.      26048304.      22253596.
 dot_product(a, b) =   7.61173658E+09
 press any key to continue

外积(三维)

program abc
  implicit none
  integer :: i
  integer, dimension(3) :: a, b
  integer, dimension(1, 3) :: c
  
  do i = 1, 3
    print '(1x, a2, i1, a2)', 'a(', i, ')='
    read *, a(i)
  end do

  do i = 1, 3
    print '(1x, a2, i1, a2)', 'b(', i, ')='
    read *, b(i)
  end do

  c = reshape((/a(2)*b(3) - a(3)*b(2), a(3)*b(1) - a(1)*b(3), a(1)*b(2) - a(2)*b(1)/), (/1, 3/))

  print *, c

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 a(1)=
3
 a(2)=
5
 a(3)=
8
 b(1)=
2
 b(2)=
2
 b(3)=
6
          14          -2          -4
 press any key to continue

矩阵

当符合矩阵各运算的定义时,有

写法意义
a + ba + b(相加)
a - ba - b(相减)
alpha*aαa(标量倍)
matmul(a, b)ab(矩阵乘法)
transpose(a)aT(转置)
a * b元素之积
a / b元素之商
a**n元素的n次方

使用reshape表示单位矩阵(3×3)

program matrix_e
  implicit none
  integer :: i
  integer, dimension(3, 3) :: p = reshape((/1, (0, i = 1, 3), 1, (0, i = 1, 3), 1/), (/3, 3/))
  !reshape(/(数或变量, 变量 = 起始值, 终止值, 递增值), (/行数, 列数/)/)
  print *, 'p = '
  print '(1x, 3i5)', p

  print*, 'press any key to continue'
  read*,

stop

end
     p =
         1    0    0
         0    1    0
         0    0    1
     press any key to continue

使用sum求行之和,列之和与总和

program jhjh
  implicit none
  integer::i,j
  real,dimension(3,3)::a
  real::sum0
  real,dimension(3)::sum1,sum2

    print*,'a11...a13='
    read*,(a(1,j),j=1,3)
    print*,'a21...a23='
    read*,(a(2,j),j=1,3)
    print*,'a31...a33='
    read*,(a(3,j),j=1,3)

    print*,'a='
    do i=1,3
      print'(1x,3f5.0)',(a(i,j),j=1,3)
    end do

    sum0=sum(a) !总和
    sum1=sum(a,dim=1) !列之和
    sum2=sum(a,dim=2) !行之和

    print'(1x,7hsum(a),f5.0)',sum0
    print'(1x,13hsum(a,dim=1),3f5.0)',sum1
    print'(1x,13hsum(a,dim=2),3f5.0)',sum2

    print*,'press any key to continue'
  read*,

stop

end
 a11...a13=
444
545
555
 a21...a23=
774
515
112
 a31...a33=
564
213
232
 a=
  444. 545. 555.
  774. 515. 112.
  564. 213. 232.
 sum(a),3954.
 sum(a,dim=1),1782.1273. 899.
 sum(a,dim=2),1544.1401.1009.
 press any key to continue

矩阵与向量的积

用到如下表现:

a(i, k: m) !矩阵a的第i行的第k到第m个元素
a(i, k: m: 2) !矩阵a的第i行的第k到第m个元素,递增量为2
a(i, :) !矩阵a的第i行
a(:, j) !矩阵a的第j列

Matrix × Column Vector

program matrix
  implicit none
  integer, dimension(3, 3) :: a
  integer, dimension(3) :: u = (/5,8,2/), v

    a = reshape((/4,5,8,7,4,1,5,6,2/), (/3, 3/))
    v = matmul(a, u)
    print '(1x, 3i6)', v

    print *,'press any key to continue'
    read *,

  stop

end
         86    69    52
     press any key to continue
Two-dimensional Rotation

program abc
  implicit none
  real :: t, x, y
  real, dimension(2, 2) :: a
  real, dimension(2, 1) :: b, c

  print *, 'x, y, t = '
  read *, x, y, t

  a = reshape((/cos(t), sin(t), -sin(t), cos(t)/),(/2, 2/))
  c = reshape((/x, y/), (/2, 1/))
  b = matmul(a, c)

  print *, 'new x, new y = ', b

  print *,'press any key to continue'
  read *,
    
  stop
    
end
 x, y, t =
222,333,0
 new x, new y =    222.000000       333.000000
 press any key to continue

Row Vector × Matrix

program matrix
  implicit none
  integer, dimension(3, 3) :: a
  integer, dimension(3) :: u = (/5,8,2/), v

    a = reshape((/4,5,8,7,4,1,5,6,2/), (/3, 3/))
    v = matmul(u, a)
    print '(1x, 3i6)', v

    print *,'press any key to continue'
    read *,

  stop

end
         76    69    77
     press any key to continue

Row Vector × Column Vector

program matrix
  implicit none
  integer :: i
  integer, dimension(3, 3) :: a
  integer, dimension(3, 1) :: u
  integer, dimension(1, 3) :: v

    u = reshape((/2, 3, 3/),(/3, 1/))
    v = reshape((/2, 3, 3/),(/1, 3/))
    a = matmul(u, v)
    print *, 'a = '
    print '(1x, 3i6)', a

    print *,'press any key to continue'
    read *,

  stop

end
     a =
          4     6     6
          6     9     9
          6     9     9
     press any key to continue

Gaussian-Jordan Elimination

Forward Elimination

For k = 1, 2, ..., n - 1 & j = k + 1, ..., n
akj = akj / akk
bk = bk / akk

For i = k + 1, ..., n & j = k + 1, ..., n
aij = aij - aikakj
bi = bi - aikbk

Back Subsitution

xn = bn / ann

For k = n - 1, n - 2, ..., 1
xk = bk - Σnj=k+1 akjbj

Example

Written Calculation

Program

ATTENTION
The program provided below can only solve full rank matrices, and it is NOT accurate enough.

program gaussian_jordan
  implicit none
  integer, parameter :: nn = 10
  integer :: i, k, kpl, n
  real :: p, q
  real, dimension(nn, nn) :: a
  real, dimension(nn) :: b, x
    !input
    print *, 'n = '
    read *, n
    print '(1x, 25hEnter_an_Augmented_Matrix)'
    do i = 1, n
      read *, a(i, 1: n), b(i)
    end do
    !forward elimination
    do k = 1, n - 1
      kpl = k + 1
      p = a(k, k)
      a(k, kpl: n) = a(k, kpl: n) / p !affecting accuracy
      b(k) = b(k) / p !affecting accuracy
      do i = kpl, n
        q=a(i, k)
        a(i, k + 1: n) = a(i, kpl: n) - q * a(k, kpl :n)
        b(i) = b(i) - q * b(k)
      end do
    end do
    !Back Subsitution
    x(n) = b(n) / a(n, n)
    do k = n - 1, 1, -1
      kpl = k + 1
      x(k) = b(k) - dot_product(a(k, kpl: n), x(kpl: n))
    end do
    !output
    print *, 'vector x = '
    print *, x(1: n)
    
    print *,'press any key to continue'
    read *,

  stop

end
     n =
    3
     Enter_an_Augmented_Matrix
    -2,1,0,-7
    5,3,2,3
    3,0,1,7
     vector x =
       1.99999988      -3.00000024       1.00000012
     press any key to continue

添加新评论

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