計算物理学II (第７回)

今回の内容• 関数副プログラム・サブルーチン副プログラム(Fortran)• 関数(C)• 何度も使う⼀連の実⾏⽂をまとめて主プログラムの外にまとめる。
• 数学の関数と考えてよいが、もっと複雑なものも関数にできる。

関数副プログラム(Fortran)• 関数：与えられた⼊⼒変数に対して操作をして値を返す⼀連の⼿続き• 例えば 4x3-5xを何度も様々なxに対して計算したいとする。
プログラムのあちこちに4x3-5xと繰り返して書くとミスのもとになる。
また後から関数形に変更が必要になったときにすべて探し出して変更する必要がある。
• 関数副プログラムとしてf(x)=4x3-5xを計算する部分を独⽴させる。
• 関数副プログラムを置く場所
1. 主プログラム(program ~end program) の中に置く(内部副プログラム)
2. 主プログラムの外(end programの後)に置く(外部副プログラム)
3. モジュールの中に置く(モジュール副プログラム)

関数副プログラム
例１：主プログラムの中に関数副プログラムを置く場合(内部副プログラム)
program example2
implicit none
double precision :: xx = 0.0d0
doprint *, x, func(x) ! ここでfunc(x)が呼び出される。関数はcontainsより下で定義されている。
x = x + 0.10d0
if( x > 5.0d0) exit
end docontains! contains以下に副プログラムを列挙する
function func(x)! 関数⾒出しdouble precision, intent(in) :: x ! xを倍精度実数として宣⾔。intent(in)がつくと関数内で値を変更できない
double precision :: func
! 関数の型を宣⾔する(倍精度実数型)
func = 4.0d0*x**3 – 5.0d0*x
return! 関数を呼び出された元のプログラムに戻る
end function func
end program example2

関数副プログラム
例2：主プログラムの外に関数副プログラムを置く場合(外部副プログラム)
program example2
implicit none
double precision :: x
double precision :: func ! 外部副プログラムは主プログラムから⾒えないので型宣⾔が必要
x = 0.0d0
doprint *, x, func(x) ! ここでfunc(x)が呼び出される。
x = x + 0.10d0
if( x > 5.0d0) exit
end doend program example2function func(x)! 関数⾒出しdouble precision, intent(in) :: x ! xを倍精度実数として宣⾔。intent(in)をつけると関数内で値を変更できない
double precision :: func
! 関数の型を宣⾔する(倍精度実数型)
func = 4.0d0*x**3 – 5.0d0*x
return! 関数を呼び出された元のプログラムに戻るend function func関数副プログラム
例3: モジュールの中に置く(モジュール副プログラム)
module mod_func! program よりも⼿前にmoduleを配置。 mod_funcは適当につけたモジュール名
implicit none
! module内でもimplicit none宣⾔
contains! containsより下にモジュール副プログラムを配置
function func(x)! 関数⾒出しdouble precision, intent(in) :: x ! xを倍精度実数として宣⾔。intent(in)をつけると関数内で値を変更できない
double precision :: func
! 関数の型を宣⾔する(倍精度実数型)
func = 4.0d0*x**3 - 5.0d0*x
return! 関数を呼び出された元のプログラムに戻る
end function funcend module mod_func
program example2
use mod_func! モジュールを使う場合はimplicit noneより⼿前でuse モジュール名 としてモジュールを参照
implicit none
double precision :: x = 0.0d0
doprint *, x, func(x) ! ここでfunc(x)が呼び出される。
x = x + 0.10d0
if( x > 5.0d0) exit
end doend program example2

関数副プログラムの例：外積を計算
program test
implicit none
double precision :: a(1:3), b(1:3), c(1:3)a(1) = 6.0d0; a(2) = 3.0d0; a(3) = 4.0d0b(1) = 3.0d0; b(2) = -2.0d0; b(3) = -4.0d0
c(1:3) = vectorproduct( a(1:3), b(1:3) ) ! 3要素の1次元配列型の関数
write(*,*) c(1:3)
contains
function vectorproduct(a, b)
double precision, intent(in) :: a(1:3), b(1:3)
double precision :: vectorproduct(1:3) ! 関数の型を要素数3の1次元配列として宣⾔
vectorproduct(1) = a(2)*b(3)-a(3)*b(2)vectorproduct(2) = a(3)*b(1)-a(1)*b(3)vectorproduct(3) = a(1)*b(2)-a(2)*b(1)returnend function vectorproductend program test

関数副プログラムの例(再帰関数)• 再帰(recurrsion)とは：関数副プログラム(やサブルーチン)の中で⾃分⾃⾝(関数副プログラムやサブルーチン)を引⽤すること。
• 再帰的に関数副プログラム(やサブルーチン)を呼び出す場合はrecursiveをつける。
• 関数の実⾏結果は関数名ではなくresult句で指定された変数に代⼊される。
フィボナッチ数列の例

recursive fuction fibonacci(n) result(fibo)
integer, intent(in) :: ninteger :: fibo! 結果を代⼊する変数
if (n<0) thenfibo = -1
returnend ifif (n == 0) thenfibo = 0else if (n == 1) thenfibo = 1
else
fibo = fibonacci(n-1) + fibonacci(n-2) ! ⾃分⾃⾝を呼び出す
end ifreturnend function fibonacciサブルーチン副プログラム• 関数副プログラムとほぼ同じ。違う点は• 値を返さない(関数には型があるがサブルーチンにはない)• ただし引数として(複数の)値を返すことができる。
• call⽂で呼び出す• サブルーチン副プログラムは1つの実⾏⽂となるが、関数副プログラムは実⾏⽂ではない• ⻑いプログラムを書く時は主プログラムに⻑々と書かずサブルーチンに分割するサブルーチン副プログラム
program example2
implicit none
double precision :: x, yx = 0.0d0
docall func(x,y)! call⽂でfuncサブルーチンを呼び出し
print *, x, yx = x + 0.10d0
if( x > 5.0d0) exit
end docontains! contains以下に副プログラムを列挙する
subroutine func(x, y)! サブルーチン⾒出し。⼊⼒も出⼒も引数に⼊れる。
implicit none
! サブルーチンではimplicit noneを使うことを推奨
double precision, intent(in) :: x
! intent(in)がつくとサブルーチン内で値を変更できない
double precision, intent(out) :: y ! intent(out)がつくとサブルーチン内で値をセットしなければならない
y = 4.0d0*x**3 – 5.0d0*x
return! 関数を呼び出された元のプログラムに戻る
end subroutine func
end program example2

サブルーチン副プログラム• サブルーチンの引数にはintent属性をつけることを推奨• intent(in) : ⼊⼒引数、サブルーチンに値が渡され、サブルーチン内で値の変更ができない• intent(out) : 出⼒変数、サブルーチンで値がセットされる変数• intent(inout) : ⼊出⼒変数、サブルーチンに値が渡され、サブルーチンで値を変更できる変数• サブルーチン内だけで使う変数には何もつけない。

サブルーチン例：外積を計算
program test
implicit none
double precision :: a(1:3), b(1:3), c(1:3)a(1) = 6.0d0; a(2) = 3.0d0; a(3) = 4.0d0b(1) = 3.0d0; b(2) = -2.0d0; b(3) = -4.0d0
call vectorproduct( a(1:3), b(1:3), c(1:3) ) ! ⼊⼒も出⼒も引数に⼊れる
write(*,*) c(1:3)
contains
subroutine vectorproduct(a, b, c)
implicit none
double precision, intent(in) :: a(1:3), b(1:3)
double precision, intent(out) :: c(1:3)c(1) = a(2)*b(3)-a(3)*b(2)c(2) = a(3)*b(1)-a(1)*b(3)c(3) = a(1)*b(2)-a(2)*b(1)returnend subroutine vectorproductend program test

変数の共有(内部副プログラムの場合)• 変数を主プログラムと副プログラム(関数やサブルーチン)で共有する• 基本は引数として渡すのが間違いが少ない
program test
implicit none
real*8 :: x0, a=1.0d0, b=2.0d0, c=3.0d0
x0 = 1.0d0
print *, x0, f(x0)! 副プログラムのxとしてx0を渡す
a = 2.0d0
print *, x0, f(x0)! aの値が変わったのでf(x0)の値が変わる
contains
function f(x)
! 2次関数を計算する関数副プログラム
real*8, intent(in) :: x ! xは引数として主プログラムから副プログラムに渡される
real*8 :: f
f = a*x**2 + b*x + c ! a, b, cは主プログラムの値を参照する
returnend function f
end program test
短いプログラムの場合は内部副プログラムが簡単に書けるが、内部副プログラムは独⽴性が低い(この副プログラムだけ他のコードにコピーしても使えない(a,b,cが中で未定義なので))

変数の共有(外部副プログラムの場合)• 基本は引数にすべて書く。引数に書かれていないものは参照できない
program test
implicit none
real*8 :: x0, a=1.0d0, b=2.0d0, c=3.0d0, f ! a,b,cの値は外部プログラムからは直接参照できない
x0 = 1.0d0
print *, x0, f(x0,a,b,c)! x0以外にa,b,cも引数として関数に渡す
a = 2.0d0
print *, x0, f(x0,a,b,c)
end program test
function f(x,a,b,c)! 様々なa,b,cの2次関数を計算をしたい場合はa,b,cも引数に⼊れる
real*8, intent(in) :: x,a,b,c
real*8 :: f
f = a*x**2 + b*x + c
returnend function f
外部プログラムは独⽴性が⾼い。計算に必要な変数がすべて引数に⼊っているので外部プログラムだけを他のコードで使うことができる。ただし引数が煩雑になる。

変数の共有(モジュール副プログラム)• モジュール内に変数を記述し、モジュールを参照する• 引数には⾊々な値で計算することを想定されるものだけを⼊れる
module func_mod !モジュールはprogramより⼿前に配置
implicit none
real*8, save :: a, b, c
!副プログラムで使う変数をsave属性つきで定義!save属性をつけておくと前に代⼊された値を保持し続ける
contains
function f(x) ! モジュール副プログラム
real*8, intent(in) :: x
real*8 :: f
f = a*x**2 + b*x + c!モジュール副プログラムはモジュール内変数を参照できる
returnend function f
end module func_mod

program test
use func_mod !implicit noneの⼿間でモジュールを参照
implicit none
real*8 :: x0
a=1.0d0; b=2.0d0; c=3.0d0 !モジュール変数に値をセットできる
x0 = 1.0d0
print *, x0, f(x0) !モジュール関数も呼び出せる
a = 2.0d0
print *, x0, f(x0)
end program test

モジュールは独⽴性が⾼い。モジュールの中で他のモジュールをuse⽂で参照することも可能
module, programの順番に書く(別のファイルに分ける場合はこの順番にコンパイル)
複数のmoduleがある場合は順番(依存関係)に注意

関数(C⾔語)• 関数：与えられた⼊⼒変数に対して操作をして値を返す⼀連の⼿続き• mainも関数。int型であれば整数値をreturnで返す。
• 例えば 4x3-5xを何度も様々なxに対して計算したいとする。
プログラムのあちこちに4x3-5xと繰り返して書くとミスのもとになる。
また後から関数形に変更が必要になったときにすべて探し出して変更する必要がある。
• 関数としてf(x)=4x3-5xを計算する部分を独⽴させる。
• 関数を置く場所
1. main関数の⼿前
2. main関数の後ろ• return⽂は関数の処理を中断してその段階の値を返す。

関数
例１：関数funcがmain関数より前に読み込まれる場合#include<stdio.h>
double func(double x){
double y;y = 4.0*x*x*x – 5.0*x;
return(y);// 倍精度実数型変数yに⼊っている値を返す}
int main(){double x;}

for(x=0.0;;){
printf(“%f %f¥n”, x, func(x)); // ここで関数にxの値を渡して呼び出す
x += 0.1;
if( x > 5.0) break;}関数
例2：関数funcがmain関数の後に読み込まれる場合#include<stdio.h>
int main(){double x;
double func(double); // 関数プロトタイプ宣⾔。mainより⼿前でもよい。
for(x=0.0;;){
printf(“%f %f¥n”, x, func(x)); // ここで関数にxの値を渡して呼び出す
x += 0.1;
if( x > 5.0) break;}}double func(double x){
double y;y = 4.0*x*x*x – 5.0*x;
return(y);// 倍精度実数型変数yに⼊っている値を返す}

再帰関数• 再帰(recursion)とは：関数で⾃分⾃⾝を引⽤すること
フィボナッチ数列

int fibonacci(int n){
int fibo;
if( n<0 ) {fibo = -1;
return(fibo);}
if( n==0 ){fibo = 0;}else if(n == 1) {fibo = 1;} else {fibo = fibonacci(n-1) + fibonacci(n-2); // ⾃分⾃⾝を呼び出す}
return(fibo);}

関数での複数の値・配列の受け渡し
関数の戻り値は⼀変数なので複数変数、あるいは配列の値を戻り値としてほしいときは引数に変数や配列のアドレスを渡し、関数の中でその変数の値を更新する。
アドレス：メモリ上での番地。変数aの場合は&aがaのアドレス。配列a[n]の場合はa、または&a[0]が先頭のアドレス
例：外積を計算#include<stdio.h>
int main(){double a[3], b[3], c[3];void vectorproduct(double [3], double [3], double [3]); // voidは戻り値がないタイプの関数
a[0] = 6.0; a[1] = 3.0; a[2] = 4.0;b[0] = 3.0; b[1] = -2.0; b[2] = -4.0;
vectorproduct( a, b, c );// 配列のアドレスを渡す。(配列名が配列のアドレス)
printf(“%f %f %f ¥n”, c[0], c[1], c[2]);return 0;}
void vectorproduct(double a[3], double b[3], double c[3]){ // double *a, double *b, double *cでもよい
c[0] = a[1]*b[2]-a[2]*b[1];c[1] = a[2]*b[0]-a[0]*b[2];c[2] = a[0]*b[1]-a[1]*b[0];return;// 関数に戻り値はない}変数の共有• 基本は引数にすべて書く。引数に書かれていないものは参照できない#include<stdio.h>
double f(double x, double a, double b, double c){return(a*x*x + b*x + c);}// 2次関数を計算する関数

int main(){double x0;x0 = 1.0;
a = 1.0; b = 2.0; c = 3.0;
printf(“%f %f¥n”, x0, f(x0,a,b,c)); // 関数にx0,a,b,cを渡す
a = 2.0;
printf(“%f %f¥n”, x0, f(x0,a,b,c)); // aの値が変わったのでf(x0,a,b,c)の値が変わる}
関数fの独⽴性が⾼い。計算に必要な変数がすべて引数に⼊っているのでこの関数だけを他のコードで使うことができる。ただし引数が煩雑になる。

変数の共有(グローバル変数)• グローバル変数をmain関数とその他の関数で共有する。
• 関数はグローバル変数に依存するため独⽴性が低くなる• 基本は引数として渡すのが間違いが少ない#include<stdio.h>
double a, b, c;
double f(double x){
return(a*x*x + b*x + c);}// グローバル変数。mainの外側で宣⾔する// 2次関数を計算する関数副プログラム// a, b, cはグローバル変数の値を参照する

int main(){double x0;x0 = 1.0;
a = 1.0; b = 2.0; c = 3.0;
printf(“%f %f¥n”, x0, f(x0)); // 関数のxとしてx0を渡す
a = 2.0;
printf(“%f %f¥n”, x0, f(x0)); // aの値が変わったのでf(x0)の値が変わる}main関数についてmain関数は整数型
正常終了は0を返し、異常終了の場合は0以外を返す
返り値はシェルから参照することができる。省略してもよい。
int main(){return (1);}// エラーの場合などは1を返す

変数$?に⼀つ⼿前のコマンドの返り値が格納される

第2回レポート• 締め切り 12⽉16⽇(⾦)• 講義資料のページに掲載(https://wwwnucl.ph.tsukuba.ac.jp/~hinohara/compphys2-22/)• 作ったプログラムと回答をpdfファイルにしてmanabaで提出。
• 講義第7回までの内容です。
• FortranまたはCで作成(C++やPythonでもOK)• 授業中・チャット・メール等での質問も歓迎します。