Fortran: 大きさ１の正方形を格子点で埋める.

Question

###前提・実現したいこと
前にもFortranのプログラムについて質問させていただきましたが、今回もお世話になります。

今写真にあるような問題をとき直しています(わからないままだったので)。

![
###発生している問題・エラーメッセージ
格子点を打ち、これに線を(gnuplotで)描かせても不完全な格子点になる.
###該当のソースコード

program lattice
  implicit none
  real(8), allocatable ::  x(:,:,:)
  real(8) :: x1 = 0.0d0, x2 = 0.0d0, x3 = 0.0d0, x4 = 0.0d0
  integer :: i,j,l, m = 2, n(2)=(/ 10, 10 /)
  real(8) dx1, dx2
  open(11, file='Ex2.15_output.d')
  allocate(x(m,n(1),n(2)))
  dx1=1.0d0/dble(n(1)-1)
  dx2=1.0d0/dble(n(2)-1)
  do j=1, n(2)
    do i=1, n(1)
      x(1,i,j)=x1
      x(2,i,j)=x2
      write(11,*) x(1,j,i), x(2,i,j)
    enddo
    write(11,*) ''
    x1=x1+dx1
    x2=x2+dx2
  enddo

  do i=1, n(1)
    do j=1, n(2)
      x(1,i,j)=x3
      x(2,i,j)=x4
      write(11,*) x(1,i,j), x(2,j, i)
    enddo
    write(11,*) ''
    x3=x3+dx1
    x4=x4+dx2
  enddo
  close(11)
end program lattice

###試したこと
いろいろ考えて、x(1,i,j)成分をi=1の場合j=1,nまで描き(点を打ち)、空行を入れる. i=2の場合j=1,nまで描き、空行を入れる...これを繰り返せばgnuplotで縦棒を描くことができるのではないかと考え、同様にしてx(2,i,j)については、iとjを入れ替えればマス目(格子点に線を描いたもの)をgnuplotに描かせることができそうとも考えました。(実際一つのdo分だけで考えてみたのですが、三角形が描かれるだけでした。)これをもとにして上のようなコードを書きました。ただ、結果としてはご覧の通りで、マス目を十分描かせることができません。

ご助言いただけるとうれしいです。以上よろしくお願いします。

###補足情報(言語/FW/ツール等のバージョンなど)
gfortran, gnuplot, 数値計算のためのFortran90/95プログラミング入門(牛島) P.63

Answer 1

座標点の生成と作図用の出力を分離すべき。
I/Oは鬼門なので分けてやる方が祟りが少ない。
ある意味、そのためにわざわざ配列を使っている。
横線だけをまず引き、次に縦線を引いてゆけば楽勝。

また座標点の生成ではdx,dyを足し合わせていっているが、これは浮動小数点数では地獄の一丁目で、避けねばならない。十進数でいえば0.333を三回足して1.0にならんと状況に陥る。どうしてもやりたいなら２のべき乗で分割しておけば避けられるが。

普通は下のように整数倍して最後に割って、浮動小数点数の誤差が溜まるのを避ける。この辺は数値計算の本の第１章あたりに必ず書いてあることなので、暇がある時に、初めの方だけ見ておくと吉。丸め誤差とか桁落ちとか、式を変形しろとグズグズ説教される。

要点を言えば、実数型は浮動小数点数で数学的実数ではないので、代数則の結合則が成り立たっていない。このため、実数の代数系で等価な式でも、コンピュータで演算すると同じ値を返さない。数学的実数の代数系で等価な式の内、コンピュータの浮動小数点数の代数系で演算して、元の式の答えに近い表現を選ぶ必要がある。

program lattice
  implicit none
  real(8), allocatable ::  x(:,:,:)
  integer :: i,j, m = 2, n(2)=(/ 10, 10 /)
  real(8) dx1, dx2
  open(11, file='Ex2.15_output.d')
  allocate(x(m,n(1),n(2)))

  do j=1, n(2)
    do i=1, n(1)
      x(1,i,j)= (1.0d0 - 0.0d0) * (i - 1) / real(n(1) - 1, kind = 8)   
      x(2,i,j)= (1.0d0 - 0.0d0) * (j - 1) / real(n(2) - 1, kind = 8)
    enddo
  enddo

  ! x-direction
  do j = 1, n(2)
    do i = 1, n(1)
      write(11, *) x(1, i, j), x(2, i, j)
    end do
    write(11, *)
  end do
  
  ! y-direction
  do i = 1, n(1)
    do j = 1, n(2)
      write(11, *) x(1, i, j), x(2, i, j)
    end do
    write(11, *)
  end do
  
  close(11)
end program lattice

plot "Ex2.15_output.d" with linespoints

Comments

[退会済みユーザー]

たびたびお世話になります。本当に助かります。

数値計算の内容については初心者なので、常識を教えていただき本当に助かります。Fortranは「いつか必要になるだろう」というモチベーションでやっていましたが、誤差についてもちゃんと考慮しておく必要があるとのご助言を頂いたので、これについても勉強していこうと思います。

I/Oについては、出来る限りシンプルに冗長にならないようにとdo文を減らしてまとめることを心がけていたのですが、ご助言通りI/Oパートは分けてかくようにします。

度々ありがとうございました。