Transcript C言語

数値解析特論
プログラミングの初歩
代数方程式の数値解法ー二分法
#include <stdio.h>
#include <math.h>
PROGRAM main
IMPLICIT NONE
void main( )
{
int n;
int n_max;
double x1;
double x2;
double f_at_x1;
double f_at_x2;
double xm;
double f_at_xm;
double x_theory;
INTEGER n
INTEGER n_max
REAL*8 x1
REAL*8 x2
REAL*8 f_at_x1
REAL*8 f_at_x2
REAL*8 xm
REAL*8 f_at_xm
REAL*8 x_theory
x_theory=1.5-sqrt(1.25);
x_theory=1.5-SQRT(1.25)
n_max=100;
x1=0.0;
x2=2.0;
f_at_x1=x1*x1-3.0*x1+1.0;
f_at_x2=x2*x2-3.0*x2+1.0;
n_max=100
x1=0.0
x2=2.0
f_at_x1=x1*x1-3.0*x1+1.0
f_at_x2=x2*x2-3.0*x2+1.0
for(n=1;n<=n_max;n++){
xm=0.5*(x1+x2);
f_at_xm=xm*xm-3.0*xm+1.0;
if(f_at_xm*f_at_x1>=0.){
x1=xm;
f_at_x1=f_at_xm;
}else{
x2=xm;
f_at_x2=f_at_xm;
}
printf("%d %e %e\n",n,xm,fabs(xm-x_theory));
}
DO n=1, n_max
xm=0.5*(x1+x2)
f_at_xm=xm*xm-3.0*xm+1.0
if(f_at_xm*f_at_x1>=0.) then
x1=xm
f_at_x1=f_at_xm
else
x2=xm
f_at_x2=f_at_xm
endif
WRITE(*,*) n, xm, ABS(xm-x_theory)
ENDDO
return;
END
}
C言語
Fortran (F90規格)
高水準言語
• 人間が容易に理解できる形態(文章，数式)で，プロ
グラムを記述するもの．
FORTRAN: 1950年代に提案．科学技術分野で多用．
間違いをおかしづらい(制限が多い)
高速な実行形式
長い歴史に伴う，過去のプログラムの蓄積
科学技術分野以外では，ほとんど使われない
代表的な規格：
FORTRAN 77 : 拡張子 (*.f)
Fortran 90 / 95 :拡張子 (*.f90) → こちらで説明します
高水準言語
• 人間が容易に理解できる形態(文章，数式)で，プロ
グラムを記述するもの．
C言語: 1972年に提案．非常に広範な分野で使用．
制限が少なく，様々な用途にプログラミング可能
(OS; UNIX を記述するために開発された経緯）
反面，間違いをおかし易い．
プログラムの構造
#include <stdio.h>
#include <math.h>
PROGRAM main
IMPLICIT NONE
void main( )
{
(省略)
(省略)
return;
}
C言語
END
Fortran (F90規格)
ヘッダーファイルの読み込み
：様々な関数が用意されているが，その定義が書かれている
ファイル（ヘッダーファイル “*.h”）を読み込む
#include <stdio.h>
#include <math.h>
→ standard I/O (Input & Output): 出力・入力に関する関数
→ 数学の関数(sin, cos, tan, exp, log,…)
プログラムの本体：main関数
＞C言語では，ソースコード中にmain関数がひとつ必要．
＞プログラムの実行時には，main関数の”{“ と”}”で囲まれた範囲が，
上から順番に実行される．
＞終了時には，”return”を書いておしまい．
プログラムの構造
#include <stdio.h>
#include <math.h>
PROGRAM main
IMPLICIT NONE
void main( )
{
(省略)
(省略)
return;
}
C言語
END PROGRAM main
Fortran (F90規格)
PROGRAM:
> プログラムの名前を定義する(任意); 上の例だと，”main”という名前にしている．
IMPLICIT NONE:
> Fortranでは，変数の種類が予約されている(暗黙の型宣言)．
ex. n, m, .. は整数(INTEGER型)と暗黙のうちに約束されている．
> 数値解析では，変数の型の混同により予期しない“計算誤差”を生じる場合がある．
> そこで，“暗黙の型宣言”を無効化して，使用する変数全てについて，
逐一定義を行うようにする（お作法です）
END:
＞プログラムの最後を教えるために必ず書く．
ソースコードのながれ
#include <stdio.h>
#include <math.h>
void main( )
{
int n;
int n_max;
double x1;
double x2;
double f_at_x1;
double f_at_x2;
double xm;
double f_at_xm;
double x_theory;
変数の型宣言
INTEGER n
INTEGER n_max
REAL*8 x1
REAL*8 x2
REAL*8 f_at_x1
REAL*8 f_at_x2
REAL*8 xm
REAL*8 f_at_xm
REAL*8 x_theory
x_theory=1.5-sqrt(1.25);
x_theory=1.5-SQRT(1.25)
n_max=100;
x1=0.0;
x2=2.0;
f_at_x1=x1*x1-3.0*x1+1.0;
f_at_x2=x2*x2-3.0*x2+1.0;
n_max=100
x1=0.0
x2=2.0
f_at_x1=x1*x1-3.0*x1+1.0
f_at_x2=x2*x2-3.0*x2+1.0
（省略）
（省略）
return;
}
PROGRAM main
IMPLICIT NONE
C言語
END
変数の型宣言
Fortran (F90規格)
>まず，以降で用いる変数について，それぞれの種類（型）を宣言する必要がある．
>変数の型宣言を行わない変数は，使用できない．
変数の型の種類
>様々な型があるが，科学技術計算の場合の基本は．．．
C言語
Fortran
整数(-1,0,1,2,…):
Int型
INTEGER型
実数(0.123など): （単精度実数）
float型
REAL型(REAL*4型)
double型
REAL*8型
（倍精度実数）
コンピュータ内での数値の表現
• 全ての情報は，メモリー（半導体）上に電気的に保存される．
>メモリー上には,”素子”が並んでいる
V1 V2
V3
V4
V5
V6
V7 V8 V9
電圧Viは高・低，ふたつの状態のみとる．
>素子の電圧の高低で情報が保存される → ２つの状態のみとる．
高電圧： “0” 低電圧: “1”
と対応（２進数）
二進数
• コンピュータ上では，10進数は2進数で表現される．
EX. 十進数”150”の二進数での表記
150 = 128 + 16 + 4 +2 = 27 + 24 + 22 + 21
= 1×27 +0×26+0×25+1×24+0×23+1×22 +1×21+0×20
2進数での表記：
1 0 0 1 0 1 1 0
V1 V2
V3
V4
V5
V6
V7 V8
(8桁の数字→8bits）
８個の素子で保存
(1桁：１素子を”bit”と呼ぶ）
Quiz: 片手の5本の指だけで，0〜31まで数えてみよう．
変数型ごとのビット数
変数の型ごとにbits数(２進数の桁)は決まっている．
桁数が長いほど，大きな数まで表現できる．
bit数
値の範囲
整数(int/INTEGER型):
32bits
-2147483648〜2147483647
実数(float/REAL型) :
32bits
0.0 と ±10-38程度〜±1038程度
実数(double/REAL*8型) :
64bits
0.0 と ±10-38程度〜±1038程度
単精度と倍精度の違い
実数の表現方法（概念的説明）：
ある実数: 3141592.65358972341414141…
指数表記: 3.14159265358972341414141…×106
指数部
仮数部
仮数部
単精度（float/REAL）：
0
0
1
指数部
．．．
0
0
．．．
1
32bits
倍精度
（double/REAL*8）：
0
0
1
仮数部と指数部を分けて，それぞれあるbits数で表現する．
仮数部
．．．
指数部
0
1
1
0
1
64bits
倍精度は仮数部の桁数が長い！
0
．．．
1
単精度と倍精度の違い
有効桁（記憶しておける仮数部の桁数）は，bits数が大きいほど長い！
有効桁数
単精度（float/REAL）：
倍精度（double/REAL*8）：
７桁
→ 3.141592
１５桁
→ 3.14159265358972
倍々精度（long double/REAL*16）： ３３桁
科学技術計算では，繰り返し計算を行う場合が多く，
有効桁数の制限による計算の誤差が蓄積しやすい．
そのため，桁数の長い “倍精度型” を使うことが暗黙の了解とされる．
誤差
演算誤差 →数値の演算による誤差
数値解析には，必ず誤差が生じる
計算誤差 →近似方法に起因する誤差
(Newton法 vs 二分法)
①丸め誤差
②打ち切り誤差
演算誤差：有効桁の存在に起因する誤差
③桁落ち誤差
④情報落ち誤差
変数の型に起因する間違い
不適切な変数型を使用することで，想定されない間違いを生じる可能性がある．
Ex1: 整数型で実数を保存
Ex2: 実数型で整数を保存
int a;
double a;
a=1.5;
a=1;
printf(“%d\n”, a);
printf(“%e\n”, a);
NG
OK
変数の型に起因する間違い
不適切な変数型を使用することで，想定されない間違いを生じる可能性がある．
Ex3: 整数型同士の計算結果が
実数に成る場合
Ex4: 整数型同士の計算結果が
実数に成る場合
int a;
int b;
int c;
int a;
int b;
double c;
a=1;
b=2;
c=a/b;
a=1;
b=2;
c=a/b;
printf(“%d\n”, c);
printf(“%e\n”, c);
NG
NG
変数の型に起因する間違い
不適切な変数型を使用することで，想定されない間違いを生じる可能性がある．
Ex5: 整数型同士の計算結果が
実数に成る場合
int a;
int b;
double c;
異なる型が演算に関わるときには，
コンピュータは型の変換を勝手に行う．
→ 予期せぬ間違いの可能性
a=1;
b=2;
c=(a/b);
printf(“%e\n”, c);
NG
変数の型に起因する間違い
異なる型を用いて計算する場合には，キャストを明示するように心がける．
キャスト: 型を変換する操作
C言語：(double)a, (int)a
Fortran：DOUBLE(a), INTEGER(a)
Ex5: 整数型同士の計算結果が
実数に成る場合
c=の計算は，実数を前提としているので，
計算に用いる整数型a,bは実数に変換して用いる．
int a;
int b;
double c;
int a;
int b;
double c;
a=1;
b=2;
c=(a/b);
a=1;
b=2;
c=((double)a/(double)b);
printf(“%e\n”, c);
printf(“%e\n”, c);
NG
OK
繰り返し計算ーforループ/DOループ
•
同じ計算を繰り返したい場合がある．
単純に書けば，
EX:
5
1
a=å
n=1 n
double a;
int n;
a=0.0;
n=1;
a=a+1./(double)n;
n=2;
a=a+1./(double)n;
n=3;
a=a+1./(double)n;
n=4;
a=a+1./(double)n;
n=5;
a=a+1./(double)n;
nを1,2,3,4,5と+1づつしながら，
同じ計算をする操作．
繰り返し計算ーforループ/DOループ
•
同じ計算を繰り返したい場合がある．
単純に書けば，
EX:
5
1
a=å
n=1 n
forループを使えば，
double a;
int n;
double a;
int n;
a=0.0;
n=1;
a=a+1./(double)n;
n=2;
a=a+1./(double)n;
n=3;
a=a+1./(double)n;
n=4;
a=a+1./(double)n;
n=5;
a=a+1./(double)n;
a=0.0;
for(n=1; n<=5; n++){
a=a+1./(double)n;
}
繰り返し計算ーforループ/DOループ
•
forループの書き方 (C言語)
for( “式A”； “条件式B”； ”式C“){
….
}
A:ループの最初に行う処理
B:ループを繰り返す条件(これが満たされる間はループは繰り返される）
C:ループ内の処理(“{…}”の部分)を１回やった後に，実行される式
for(n=1； n<=5； n++){
….
}
A:最初にn=1を代入．
B: nは5以下(“<=“は”以下”という意味)の間繰り返す
C: nを1づつ増やす(“n++”は”n=n+1”という意味)
繰り返し計算ーforループ/DOループ
•
DOループの書き方 (FORTRAN)
DO n=1, 5
….
END DO
FORTRANでは，n=1,2,3,4.. のように１ループごとに，nは1づつ増える。
条件分岐
•
計算された変数の値によって，異なる処理を実行したい場合がある→条件分岐
「if 文」
FORTRAN
C言語
if(“条件式”) {
“処理A”
}else{
“処理B”
}
if(“条件式”) then
“処理A”
else
“処理B”
endif
“条件式”が成り立つ場合 → “処理A”が実行
それ以外の場合
→ “処理B”が実行
条件分岐
•
else文の部分は省略可能
「if 文」
C言語
if(“条件式”) {
FORTRAN
if(“条件式”)
then
“処理A”
“処理A”
}
endif
“条件式”が成り立つ場合 → “処理A”が実行
それ以外の場合
→ なにもしない
条件式の種類
•
C言語の場合
条件式
aとbが等しければ
aとbが等しくなければ
if(a==b) {
if(a!=b) {
aがbより大きければ
aがb以上であれば
if(a>b) {
if(a>=b) {
aがbより小さければ
aがb以下であれば
if(a<b) {
if(a<=b) {
二つ以上の条件の組み合わせ
“条件A”と”条件B”が共に成り立つ(and)ならば, if( “条件A” && “条件B”){
例: if((a==b)&&(c<=d)) {
aがbと等しく，cはd以下ならば，
“条件A”と”条件B”の何れかが成り立つ(or)ならば, if( “条件A” || “条件B”){
例： if((a==b)||(c<=d)){
aがbと等しいか，cはd以下の何れかであれば，
条件式の種類
•
FORTRANの場合
条件式
aとbが等しければ
aとbが等しくなければ
if( a==b ) then 或は if( a .eq. b) then
if(a/=b) then 或は if( a .ne. b) then
aがbより大きければ
aがb以上であれば
if(a>b) then
或は if( a .gt. b) then
if(a>=b) then 或は if( a .ge. b) then
aがbより小さければ
aがb以下であれば
if(a<b) then
或は if( a .lt. b) then
if(a<=b) then 或は if( a .le. b) then
二つ以上の条件の組み合わせ
“条件A”と”条件B”が共に成り立つ(and)ならば, if( “条件A” .and. “条件B”)
例: if((a==b) .and. (c<=d)) then
aがbと等しく，cはd以下ならば，
“条件A”と”条件B”の何れかが成り立つ(or)ならば, if( “条件A” .or. “条件B”)
例： if((a==b) .or. (c<=d)){
aがbと等しいか，cはd以下の何れかであれば，
配列ー変数の型の拡張
• 複数の変数を使用したい場合
• for/DOループにより，複数データについて同一処理をしたい
場合に有効
例：１次元配列
A1,A2,A3,．．．,A100と，100個の変数を使用し,
A1=1./1., A2=1./2., A3=1./3., …, A100=1./100.
の計算を行う場合．
double A1;
double A2;
…
double A100;
配列を使用しない場合；
配列ー変数の型の拡張
• １次元配列ーC言語の場合
配列を使用しない場合
配列を使用した場合
…
double A1;
double A2;
…
double A100;
…
double A[100];
A1=1./1.;
A2=1./2.;
…
A99=1./99.;
A100=1./100.
A[0]=1./1.;
A[1]=1./2.;
…
A[98]=1./99.;
A[99]=1./100.
配列ー変数の型の拡張
• １次元配列ーFORTRANの場合
配列を使用しない場合
配列を使用した場合
…
REAL*8 A1
REAL*8 A2
…
REAL*8 A100
…
REAL*8 A(100)
A1=1./1.
A2=1./2.
…
A99=1./99.
A100=1./100.
A(1)=1./1.
A(2)=1./2.
…
A(99)=1./99.
A(100)=1./100.
配列ー変数の型の拡張
• １次元配列
変数の型宣言
0〜99個の100個の変数
C言語： double A[100]
int N[50]
FORTRAN: REAL*8 A(100)
INTEGER N(50)
A[0], A[1], A[2], … , A[98], A[99]
N[0], N[1], N[2], … , N[48], N[49]
1〜100個の100個の変数
A(1), A(2), A(3), … , A(99), A(100)
N(1), N(2), N(3), … , N(49), N(50)
配列ー変数の型の拡張
• 配列は for/DO ループ と組み合わせると効果絶大
配列を使用しない場合
配列を使用した場合
…
double A1;
double A2;
…
double A100;
…
double A[100];
A1=1./1.;
A2=1./2.;
…
A99=1./99.;
A100=1./100.
A[0]=1./1.;
A[1]=1./2.;
…
A[98]=1./99.;
A[99]=1./100.
配列を使用した場合+for文
int i;
double A[100];
for(i=0;i<=99;i++){
A[ i ] = 1./(double)i;
}
配列の番号として，整数の変数
を代入することができる．
配列ー変数の型の拡張
• 配列は for/DO ループ と組み合わせると効果絶大
配列を使用しない場合
配列を使用した場合
…
REAL*8 A1
REAL*8 A2
…
REAL*8 A100
…
REAL*8 A(100)
A1=1./1.
A2=1./2.
…
A99=1./99.
A100=1./100.
A(1)=1./1.
A(2)=1./2.
…
A(99)=1./99.
A(100)=1./100.
配列を使用した場合+for文
INTEGER i
REAL*8 A(100);
DO i = 1, 100
A ( i ) = 1./DOUBLE(i);
END DO
配列ー変数の型の拡張
• 多次元配列
2次元配列
C： double A[10][30];
A[0][0], A[0][1], A[0][2], …, A[0][29]
A[1][0], A[1][1], A[1][2], …, A[1][29]
A[i][j]: i=0,1,2…,9
j=0,1,2,…,29
の10×30個の変数
A[2][0], A[2][1], A[2][2], …, A[2][29]
….
A[9][0], A[9][1], A[9][2], …, A[9][29]
FORTRAN： REAL*8 A(10, 30)
A ( i, j ): i=1,2,3,…,10
の10×30個の変数
j=1,2,3,…,30
A(1,1), A(1,2), A(1,3), …, A(1,30)
A(2,1), A(2,2), A(2,3), …, A(2,30)
A(3,1), A(3,2), A(3,3), …, A(3,30)
….
A(10,1), A(10,2), A(10,3), …, A(10,30)
配列ー変数の型の拡張
• 多次元配列
3次元配列
C： double A[10][30][40];
FORTRAN： REAL*8 A(10, 30, 40)
4次元配列
C： double A[10][30][40][5];
FORTRAN： REAL*8 A(10, 30, 40, 5)
．．．たしか，７次元までOKだったと思う。
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
C言語： printf関数を用いる．
int n;
double a;
double b;
n=1;
a=0.;
b=2.;
printf( “ %d %e %e \n” , n, a, b);
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
C言語： printf関数を用いる．
printf( “ %d %e %e \n” , n, a, b);
画面上への出力の様子
書式指定子
%d ：ここに整数を出力
%e ：ここに実数(指数表示; 1.0000e+3)を出力
%f ：ここに実数(小数点表示; 1000.0000)を出力
%s ：ここに文字を出力
\n :ここで改行
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
C言語： printf関数を用いる．
ここでの並びの順番で”…”内の
%d or %e の位置に値が出力．
n=1; a=0.; b=2.;
printf( “ %d %e %e \n” , n, a, b);
実際の画面上では，
1
0.0000E+00
2.00000E+00
と出力される．
書式指定子の個数と，変数の数は一致しなければならない
書式指定子の型と，変数の型は一致しなければならない
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
C言語： printf関数を用いる．
n=1; a=0.; b=2.;
printf( “ %d %e\n %e \n” , n, a, b);
実際の画面上では，
1 0.0000E+00
2.00000E+00
と出力される．
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
C言語： printf関数を用いる．
n=1; a=0.; b=2.;
printf( “v_n=%d v_a=%e v_b=%e \n” , n, a, b);
実際の画面上では，
v_n=1 v_a=0.0000E+00 v_b=2.00000E+00
と出力される．
出力ー画面（コンソールorターミナル上へ）
• 数値の出力
FORTRAN： WRITE関数を用いる．
INTEGER n
REAL*8 a
REAL*8 b
n=1
a=0.0
b=2.0
WRITE( * , * ) n, a, b
ファイルへの出力
• C言語，FORTRANともに
① ファイルを開く
② 書き出す(fprintf関数 or WRITE関数)
③ ファイルを閉じる
の３段階で行う．
出力ーファイルへ
• C言語
int n;
double a;
double b;
FILE *fp;
n=1; a=0.; b=2.;
fp=fopen ( “output.txt”, “w”);
fprintf(fp, “ %d %e %e \n” , n, a, b);
fclose(fp);
出力ーファイルへ
• C言語
FILE *fp;
ファイルを区別するための変数（お作法）
fp = fopen ( “ output.txt ”, “w”);
output.txtという名前のファイルを，書出しのために開く．
“w” : 書き出し
output.txtというファイルが無ければ，新たに作り，
ファイルが有れば，それを消去した上で，新規に作成
“a” : 追加書出し
既存のファイル(output.txt)の末尾に追加して出力する．
出力ーファイルへ
• C言語
fprintf ( fp, “ %d %e %e \n” , n, a, b);
実際の書出しには，printfの代わりに，fprintf 関数
第一引数に，ファイル変数”fp”を指定する．
fclose (fp); ファイルを閉じる（出力終了）
出力ーファイルへ
• FORTRAN
INTEGER n
REAL*8 a
REAL*8 b
(省略)
OPEN(UNIT=10, FILE=’output.txt', STATUS='replace')
WRITE(10, *) n, a, b
CLOSE( 10 )
出力ーファイルへ
• FORTRAN
OPEN(UNIT=10, FILE=’output.txt', STATUS='replace’)
装置番号でファイルを区別する（UNIT=10 なら，装置番号１０）
STATUS=xxx は，開くファイルの状態
‘replace’ : 既にあるファイルを置き換える．無ければ，新規に
作成する．
‘new’ : 新規に作成する（無ければエラー）
‘old’ : 既にあるファイルを読み込みなどのために開く
出力ーファイルへ
• FORTRAN
OPEN(UNIT=10,
WRITE(10,
*) n, FILE=’output.txt',
a, b
STATUS='replace’)
装置番号１０番のファイルへ出力
OPEN(UNIT=10, FILE=’output.txt', STATUS='replace’)
CLOSE(10)
装置番号１０番のファイルを閉じる．