数値解析法及び演習 第三回

DOループの応用
  1. 今日の目標
  2. 異なる型の混合演算について
  3. gnuplotの基本
  4. 配列
  5. 最小二乗法
  6. 今日の宿題

[1] 今日の目標

  1. 数の型の理解を深める
  2. 配列の概念を学ぶ

[2] 異なる型の混合演算について

FORTRANでは実数と整数を区別して変数に格納します. 正しい型の変数を用意しておかないと計算結果がおかしくなります. 変数の型には注意をはらうこと!

  • 実数は必ず小数点がつく! 3.D0は実数，3は整数． 整数型と実数型の変数を混ぜて使うときには注意が必要である．以下に考えられるパ ターンを挙げる． 
          program example1
      implicit none
      integer :: m, n, l
      double precision :: a, b, c, d, s
      real :: f, r

      l=1/2
      a=1/2
      
      m=1.D0/2.D0
      b=1.D0/2.D0
      
      n=5/2.D0
      c=5/2.D0
      
      r = 16.1
      s = 3.D0

      s = s + r
      
      write(*,'(3(I2,x),2(E24.16,x))') l, m, n, a, b
      write(*,'(2(E24.16,x))') c, d
      write(*,'(2(E24.16,x))') s, r 

      end program example1
    結果はどうなったであろうか? dの値は？

    • [3] gnuplotの基本

    gnuplot(ニュープロットとよむ)は，UNIXシステムにおける簡便なグラフ描画 ソフトである． 起動は 
    gnuplot
    で行われる．
    最も基本的なデータは 
    -1.0   -3.2
 0.0   -0.9
 1.0    1.2
 2.0    2.8
 3.0    5.0
    といった形で与えられる．左側のデータをx座標，右側のデータ(ファイル名を plot.datとする)をy座標とおもってxy平面上にプロットするときは，gnuplotを起動した画面で 
    plot 'plot.dat'
    とするとグラフが描かれる．さらに 
    plot 'plot.dat' w l
    とするとデータ点を線で結ぶ．また，gnuplotにはいろいろな関数が用意されていて， 
    plot sin(x)
    とするとサインカーブを描いてくれる．三次元プロットも可能で 
    sp x**2 + y**2
    とすると　z=x**2 + y**2を描画する．

    [4] 配列

    通常の変数（整数型，実数型ともに）は数値を一つしか格納できない．これでは多くのデータを扱うときに困る．そこで「配列」とよばれるものが用意されている．配列は複数の要素からなり，要素の番号を指定することで複数の数値を格納することができる．

    配列名は英字で始まり，31文字以内の任意の英数字およびアンダースコア（＿）からなる．大文字と小文字は区別されない．

    配列の宣言:
    宣言文 配列名（寸法の下限：寸法の上限）

    例1 
    double precision :: a(10)
    a(1)からa(10)まで要素を10個持つ倍精度実数型配列を宣言

    例2 
    integer :: i(-1:8)
    i(-1)からi(8)まで要素を10個もつ整数型配列を宣言

    例： 
           program example2
       implicit none 
       integer :: i
       double precision :: a(10)

       a(1)=1.D0
       a(2)=2.D0
       a(3)=3.D0


       do i=1,10
          write(*,*) i,a(i)
       enddo

       end program example2
    a(4)以降の要素はどう表示されたでしょうか？

    練習問題1

    次のプログラムを改造して，

    それぞれの総和（計3つ）を求めなさい．なお，以下のプログラムを基にすること．それぞれの総和を，sum(1), sum(2), sum(3)に格納すること． 

           program example3
       implicit none 
       integer :: i
       double precision :: a(10),b(10),sum(3)

       do i=1,3
          sum(i)=0.D0  ! 三つの和をゼロにリセットする
       enddo
    
       do i=1,10
          a(i)=1.D0*i  ! a(i)には実数の1.D0, 2.D0,... 10.D0が格納される
          b(i)=2.D0*i  ! b(i)には実数の2.D0, 4.D0,... 20.D0が格納される
       enddo

       **ここに必要なプログラムを書き足して下さい**

       end program example3

    練習問題2

    上のプログラムを改造して，

    を求めなさい．

    [5] 最小二乗法

    二つの量を測定したときに，そのあいだに線形の関係が期待されることがよくある． もっともありそうな傾き(A)と切片(B)を求めるのが最小二乗法である．

    最小二乗法では, 測定されたデータの組(x,y)にフィットする直線を求める. データxは正確に求められているとし, またデータyのバラツキは正規分布をしているものとする. この場合, フィットする直線をy=Ax+Bとすると, 傾きA,切片Bはそれぞれ（ここでΣは1からNまでの和）

     
    $\displaystyle B$ $\textstyle =$ $\displaystyle {N\sum x_iy_i-\sum x_i\sum y_i\over N\sum x_i^2-(\sum x_i)^2}$  
    $\displaystyle A$ $\textstyle =$ $\displaystyle {\sum x_i^2\sum y_i-\sum x_i\sum x_i y_i\over N\sum x_i^2-(\sum x_i)^2}$

    となる. したがって切片Aと傾きBを算出するためには次の和をもとめればよい．

    配列へのデータの読み込み

    先ほどgnuplotで表示させたデータ plot.dat を配列x(5),y(5)に読み込ませる．read文とdoループを組み合わせて以下のようにかける．以下のプログラムを実行し，正しく配列にデータを読み込んでいることを確認すること．

           implicit none 
       integer :: i
       double precision :: x(5), y(5)
       
       open(1,file='plot.dat')
       
       do i=1,5
         read(1,*) x(i), y(i)
         write(*,*) x(i), y(i)
       enddo
       end

    練習問題3

    先ほどのプログラムを改造し，最小二乗法によりフィットした直線の傾きAと切片Bを求めるプログラムを作成せよ. AとBを求めたら, その直線とデータをgnuplotで描き, フィッテイングが正しいことを確かめよ. データと直線を重ねてプロットするには 
    plot 'plot.dat', 3.4 * x + 2.2
    のように，コンマで続ける．2.2や3.4は自分で求めたA, Bの値に変更すること．

    [6] 今日の宿題

    練習問題3のプログラム, 答えを送ること. 〆切は10/25.


    日程表へ戻る <<