”’lŒvŽZ‚Ì‚½‚߂̃ƒ[ƒbƒ^ƒXƒg[ƒ“i–¢Š®j

Numerical Rosetta Stone

1

目的

１．数種のプログラム言語を用いて，数値計算プログラムの例を提供する． ２．これによって新しい言語を学習する際に，既知の言語の知識が利用できる． ３．言語はC, Java, Fortran, Perl, Ruby, Python の６つを取り上げた．

４．前３つはコンパイラ言語，後３つはスクリプト言語である．好みと状況に応じて使 い分けるのが良いだろう．

[例題１]（Hello world）

コマンドラインに文字列を出力する定番プログラム．

C:

/* hello.c 実行は $> gcc hello.c の後 $> ./a.out とする */ #include <stdio.h>

main(){

printf("Hello world of C!\n"); }

Java:

// hello.java 実行は $> javac hello.java の後 $> java hello とする public class hello{

public static void main(String args[]){

System.out.println("Hello world of Java!\n"); }

}

Perl:

# hello.pl 実行は $> perl hello.pl とする print "Hello world of Perl!\n";

Ruby:

# hello.rb 実行は $> ruby hello.rb とする print "Hello world of Ruby!\n"

Python

# hello.py 実行は $> python hello.py とする print ’Hello world of Python!’

[類題１] 円周率

言語系既存の機能を用いて円周率の値を入手する．これを使えばπ = 3.141592· · · とい

う数値を覚えなくても済む．

C:

/* pi.c 実行は $> gcc pi.c -lm の後 $> ./a.out とする */ #include <stdio.h> #include <math.h> main(){ double pi=4*atan(1); printf("%lf\n", pi); # あるいは pi=M_PI; printf("%lf\n", pi); # としてもよい } Java:

// pi.java 実行は $> javac pi.java の後 $> java pi とする public class pi{

public static void main(String args[]){ double pi=4*Math.atan(1);

System.out.println(pi); }

} Perl:

# pi.pl 実行は $> perl pi.pl とする $pi=4*atan2(1,1);

print $pi; Ruby:

# pi.rb 実行は $> ruby pi.rb とする print 4*Math.atan(1) # あるいは include Math print PI # としても良い（PI と E は予約定数） Python:

2

逐次計算

[例題２] フィボナッチ数列 漸化式 Fn+1 = Fn+ Fn−1 を初期条件 F0 = 1, F1 = 1 で解いた数列をフィボナッチ数 列と呼ぶ．この数列は数学・物理の色々な場面で登場する． C: /* fibo.c */ #include <stdio.h> int fibo(int n){ int j, f0=1, f1=1, ff; if(n==0) ff=f0; elseif(n==1) ff=f1; else{ for(j=1; j<n; j++){ ff=f0+f1; f0=f1; f1=ff; } ff=f1; } return ff; } main(){ int j, ans, n=50; for(j=0; j<=n; j++){ ans=fibo(j); printf("F(%d)=%d\n", j, ans); } } Java: // fibo.java

public class fibo{

public static int F(int n){ int j, f0=1, f1=1, ff=2; if(n==0) ff=f0; else if(n==1) ff=f1; else{ for(j=1; j<n; j++){ ff=f0+f1; f0=f1; f1=ff; } }

return ff; }

public static void main(String args[]){ int j, ans, n=50; for(j=0; j<=n; j++){ ans=F(j); System.out.println("F("+j+")="+ans); } } } Perl: # fibo.pl $n=50; for($j=0; $j<=$n; $j++){ $ff=fibo($j); print "F($j)=$ff \n"; } sub fibo{ $f0=1; $f1=1; if($_[0]==0){ $ff=$f0; } elsif($_[0]==1){ $ff=$f1; } else{ for($k=1; $k<$_[0]; $k++){ $ff=$f0+$f1; $f0=$f1; $f1=$ff; } } return $ff; } Ruby: #fibo.rb def fibo(n) f0, f1=1, 1 if n==0 then ff=f0 elsif n==1 then ff=f1 else

for j in 1..(n-1) ff=f0+f1 f0, f1=f1, ff end end return ff end for k in 0..50 print("F(",k,")=",fibo(k),"\n") end Python: # fibo.py def fibo(n): f0, f1=1, 1 if n==0: ff=f0 elif n==1: ff=f1 else: for j in range(1, n): ff=f0+f1 f0, f1=f1, ff return ff for k in range(50): print ’F(’,k,’)=’,fibo(k) （注）得られた結果が本当に正しいかどうかにも注意せよ．間違っている場合には、その 理由を考えてみよ． またサブルーチンを再帰的な記述 fibo(n)=fibo(n-1)+fibo(n-2) に直して，実行速度を較べてみるとおもしろい．次の類題「階乗 n! の計算」はそういう 再帰アルゴリズムの典型として有名である．

[類題２] 階乗n! の計算 C: /* factorial.c */ #include <stdio.h> int factorial(int n){ int ff; if(n==0) ff=1; else{ ff=(n*factorial(n-1)); } return ff; } main(){ int j; for(j=0; j<=20; j++){ printf("%d!=%d\n", j, factorial(j)); } } Java: // factorial.java

public class factorial{

public static int F(int n){ int j, ff; if(n==0) ff=1; else{ ff=n*F(n-1); } return ff; }

public static void main(String args[]){ int j, ans, n=20; for(j=0; j<=n; j++){ ans=F(j); System.out.println(j+"!="+ans); } } } Perl: # factorial.pl

$n=20; for($j=0; $j<=$n; $j++){ $ff=factorial($j); print "$j!=$ff \n"; } sub factorial{ if($_[0]==0){ $ff=1; } else{ $ff=$_[0]*factorial($_[0]-1); } return $ff; } Ruby: #factorial.rb def factorial(n) if n==0 then ff=1 else ff=n*factorial(n-1) end return ff end for j in 0..20 print(j,"!=",factorial(j), "\n") end Python def factorial(n): if n==0: ff=1 else: ff=n*factorial(n-1) return ff for j in range(20): print j,’!=’,factorial(j)

3

求解法

[例題３] ニュートン法 方程式 f (x) = 0を満たす変数x を求める求解法のひとつ，ニュートン法は，適当な初 期値 x0 から始めて，漸化式 xn+1= xn− f (xn) f0(xn) の繰り返しにより収束値を計算するという方法である．n→ ∞ で収束する場合，上式は f (x_∞) = 0 に帰着するから，x_∞ は f (x)のゼロ点となる． 関数 f (x) = x2− 2として，２の平方根√2 = 1.41421· · · を求めてみよう． C: /* newton.c */ #include <stdio.h>

double rhs(double x, double a){ return ((x+a/x)/2); } main(){ int n; double x, a; a=2; x=a; for(n=0; n<10; n++){ x=rhs(x, a); printf("%lf\n", x); } } Java: // newton.java

public class newton{

public static double rhs(double x, double a){ return ((x+a/x)/2);

}

public static void main(String args[]){ int n; double x, a; a=2; x=a; for(n=0; n<10; n++){ x=rhs(x,a); System.out.println(x);

} } } Perl: # newton.pl $a=2; $x=$a; for($n=0; $n<10; $n++){ $x=rhs($x, $a); print "$x \n"; } sub rhs{ return (($_[0]+$_[1]/$_[0])/2); } Ruby: #newton.rb def rhs(x, a) return (x+a/x)/2 end a=2 x=a for n in 1..10 x=rhs(x, a) p x end Python: #newton.py def rhs(x, a): return ((x+a/x)/2) a=2.0 x=a for n in range(10): x=rhs(x, a) print x

[類題３] ２分法 もうひとつの求解法である ２分法 も取り上げよう．こちらは，連続関数 f (x) に対す る中間値の定理 x0< x1 のとき f (x0)· f(x1) < 0 ならば， x0 < x < x1 でf (x) = 0 を満たすx が少なくとも一つ存在する を利用して，中点における関数の正負を調べて範囲を狭めていくという方法である． やはり関数f (x) = x2− 2として，平方根 √2 を計算してみよう． C: /* bisection.c */ #include <stdio.h>

double residue(double x, double a){ return (x*x-a);

}

main(){ int n;

double x, a, left, right, center; a=2.0; left=0; right=a; center=(left+right)/2; for(n=0; n<30; n++){ if(residue(center, a)>0){ right=center; }else{ left=center; } center=(left+right)/2; printf("%15.13f\n", center); } } Java: // bisection.java

public class bisection{

public static double residue(double x, double a){ return (x*x-a);

}

public static void main(String args[]){ int n;

double x, a, left, right, center; a=2.0;

left=0; right=a; center=(left+right)/2; for(n=0; n<30; n++){ if(residue(center, a)>0){ right=center; }else{ left=center; } center=(left+right)/2; System.out.println(center); } } } Perl: # bisection.pl $a=2.0; $left=0.0; $right=$a; $center=($left+$right)/2.0; for($n=0; $n<30; $n++){ if(residue($center, $a)>0){ $right=$center; }else{ $left=$center; } $center=($left+$right)/2.0; print "$center\n"; } sub residue{ return ($_[0]*$_[0]-$_[1]); } Ruby: #bisection.rb def residue(x, a) return (x*x-a) end a=2.0 left, right=0, a center=(left+right)/2; for n in 1..30

if residue(center, a)>0 then right=center else left=center end center=(left+right)/2 p center end Python: #bisection.py def residue(x, a):

return (x*x-a)

a=2.0 # What happen if you write "=2" ? left, right=0, a center=(left+right)/2 for n in range(30): if residue(center, a)>0: right=center else: left=center center=(left+right)/2 print center

4

定積分

[例題４] 定積分 有限区間の定積分を台形公式で数値計算する． ∫ b a f (x)dx = N∑−1 j=0 h 2 · (f(xj) + f (xj+1)) , ( h = b− a N , xj = a + j· h ) 0≤ x ≤ 1, f(x) = 1/(x2+ 1)の場合，結果は π/4 となるはずである． C: /* sekibun.c */ #include <stdio.h> #include <math.h> double f(double x){ return (1/(x*x+1)); } main(){ int j, N=100; double h, x, S; h=1/(double)N; S=0.0; for(j=0; j<N; j++){ x=h*(double)j; S+=(h/2)*(f(x)+f(x+h)); } printf("%lf %lf\n", S, atan(1)); } Java: // sekibun.java

public class sekibun{

public static double f(double x){ return (1/(x*x+1));

}

public static void main(String args[]){ int j, N; double h, x, S; N=100; h=1/(double)N; S=0.0; for(j=0; j<N; j++){

S+=(h/2)*(f(x)+f(x+h)); } System.out.println(S); System.out.println(Math.atan(1)); } } Perl: # sekibun.pl $N=100; $h=1.0/$N; $S=0.0; for($j=0; $j<$N; $j++){ $x=$h*$j; $S+=($h/2)*(f($x)+f($x+$h)); } print "$S\n"; print atan2(1,1); sub f{ return (1/($_[0]*$_[0]+1)); } Ruby: # sekibun.rb def f(x) return (1/(x*x+1)) end N=100 h=1.0/N s=0.0 # 変数には小文字を使う！ for j in 1..N x=h*(j-1) s+=(h/2)*(f(x)+f(x+h) end p s p Math.atan(1.0)) Python: # sekibun.py import math def f(x):

return (1/(x*x+1)) N=100 h=1.0/N S=0.0 for j in range(N): x=h*j S+=(h/2)*(f(x)+f(x+h)) print S print math.atan(1) [類題４] 完全楕円積分 第１種完全楕円積分は K(k) = ∫ π 2 0 dx √ 1− k2_sin2_x (0≤ k ≤ 1) で定義される．kの関数としてK(k)の数値を得るには，ループを２重にすれば良い．縦 軸を K(k)，横軸を kとしてグラフを描いてみよ． C: #include <stdio.h> #include <math.h>

double f(double x, double k){

return (1/(sqrt(1-pow(k*sin(x),2)))); } main(){ int j, n, N=100; double x, k; double pi, h, S; k=0.0; pi=4*atan(1.0); h=pi/(2*(double)N); for(j=0; j<50; j++){ S=0.0; for(n=0; n<N; n++){ x=h*(double)n; S+=(h/2)*(f(x, k)+f(x+h, k)); } printf("%lf %lf\n", k, S); k+=0.02; }

} Java:

// elliptic.java

public class elliptic{

public static double f(double x, double k){

return (1/(Math.sqrt(1-Math.pow(k*Math.sin(x),2)))); }

public static void main(String args[]){ int i, j, N;

double pi, h, x, k, S; N=100;

pi=4*Math.atan(1); h=(pi/2)/(double)N; for(i=0; i<50; i++){

k=0.02*(double)i; S=0.0; for(j=0; j<N; j++){ x=h*(double)j; S+=(h/2)*(f(x,k)+f(x+h,k)); } System.out.println(k+" "+S); } } } Perl: # elliptic.pl $N=100; $pi=4*atan2(1,1); $h=($pi/2)/$N;

for($i=0; $i<50; $i++){ $k=0.02*$i; $S=0.0; for($j=0; $j<$N; $j++){ $x=$h*$j; $S+=($h/2)*(f($x,$k)+f($x+$h,$k)); } print "$k $S\n"; } sub f{

return (1/(sqrt(1-($_[1]*sin($_[0]))**2))); } Ruby: #elliptic.rb def f(x, k) return (1/Math.sqrt(1-(k*Math.sin(x))**2)) end N=100 pi=4.0*Math.atan(1.0) h=(pi/2)/N for i in 0..50 k=i/50.0 s=0.0 for j in 1..n x=h*(j-1) s+=(h/2)*(f(x, k)+f(x+h, k)) end print(k, " ", s, "\n") end Python: # elliptic.py import math def f(x, k): return (1/(math.sqrt(1-pow(k*math.sin(x),2)))) N=100 pi=4*math.atan(1) h=(pi/2)/N for i in range(50): k=0.02*i S=0.0 for j in range(N): x=h*j S+=(h/2)*(f(x,k)+f(x+h,k)) print k,’ ’,S

5

常微分方程式

[例題５] ロレンツ方程式 カオス現象の発見につながったロレンツ方程式は dx dt =−σ(x − y) σ = 10 dy dt =−y + rx − xz r = 28 dz dt =−bz + xy b = 8/3 で与えられる．適当な初期条件の下でこれらを解き，グラフを描け． 以下にはルンゲ・クッタ法の場合を示すが，他にも修正オイラー法やホイン法などがあ る．これらは，連立微分方程式 dx dt = f (x, t) を t = nhと差分化するとき ホイン法： k1 = h· f(xn, t), k2= h· f(xn+ k1, t + h) ⇒ xn+1= xn+ 1 2(k1+ k2) ルンゲ・クッタ法： k1 = h· f(xn, t), k2 = h· f(xn+ 1 2k1, t + h 2), k3 = h· f(xn+ 1 2k2, t + h 2), k4 = h· f(xn+ k3, t + h) ⇒ xn+1= xn+ 1 6(k1+ 2· k2+ 2· k3+ k4) と，陽解法（explicit method）で計算される． C: /* lorenz.c */ #include <stdio.h> #define sigma 10.0 #define r 28.0 #define b 8.0/3.0

double X(double x, double y, double z){ return (-sigma*(x-y));

}

return (-y+r*x-x*z); }

double Z(double x, double y, double z){ return (-b*z+x*y);

}

void rk4(double x,double y,double z, double *xx,double *yy,double *zz,double h){ double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; x1=h*X(x,y,z); y1=h*Y(x,y,z); z1=h*Z(x,y,z); x2=h*X(x+x1/2,y+y1/2,z+z1/2); y2=h*Y(x+x1/2,y+y1/2,z+z1/2); z2=h*Z(x+x1/2,y+y1/2,z+z1/2); x3=h*X(x+x2/2,y+y2/2,z+z2/2); y3=h*Y(x+x2/2,y+y2/2,z+z2/2); z3=h*Z(x+x2/2,y+y2/2,z+z2/2); x4=h*X(x+x3,y+y3,z+z3); y4=h*Y(x+x3,y+y3,z+z3); z4=h*Z(x+x3,y+y3,z+z3); *xx=x+(x1+2*x2+2*x3+x4)/6; *yy=y+(y1+2*y2+2*y3+y4)/6; *zz=z+(z1+2*z2+2*z3+z4)/6; } main(){ int j, N=2000; double h, x, y, z, xx, yy, zz; /* initial condition */ x=1.0; y=1.0; z=1.0; printf("%lf %lf %lf\n", x, y, z); /* time step */ h=0.01; /* time development */ for(j=0; j<N; j++){ rk4(x,y,z,&xx,&yy,&zz,h); x=xx; y=yy; z=zz; printf("%lf %lf %lf\n", x, y, z); } } Java:

// lorenz.java

public class lorenz{

static double sigma=10.0, r=28.0, b=8.0/3.0;

public static double X(double x,double y,double z){ return (-sigma*(x-y));

}

public static double Y(double x,double y,double z){ return (-y+r*x-x*z);

}

public static double Z(double x,double y,double z){ return (-b*z+x*y);

}

public static void main(String args[]){ int j, N=2000; double h, x, y, z; double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4; // initial condition x=1.0; y=1.0; z=1.0; System.out.println(x+" "+y+" "+z); // time step h=0.01; // time development for(j=0; j<N; j++){ x1=h*X(x,y,z); y1=h*Y(x,y,z); z1=h*Z(x,y,z); x2=h*X(x+x1/2,y+y1/2,z+z1/2); y2=h*Y(x+x1/2,y+y1/2,z+z1/2); z2=h*Z(x+x1/2,y+y1/2,z+z1/2); x3=h*X(x+x2/2,y+y2/2,z+z2/2); y3=h*Y(x+x2/2,y+y2/2,z+z2/2); z3=h*Z(x+x2/2,y+y2/2,z+z2/2); x4=h*X(x+x3,y+y3,z+z3); y4=h*Y(x+x3,y+y3,z+z3); z4=h*Z(x+x3,y+y3,z+z3); x+=(x1+2*x2+2*x3+x4)/6; y+=(y1+2*y2+2*y3+y4)/6; z+=(z1+2*z2+2*z3+z4)/6; System.out.println(x+" "+y+" "+z); } }

} Perl: # lorenz.pl $sigma=10.0; $r=28.0; $b=8.0/3.0; sub X{ return (-$sigma*($_[0]-$_[1])); } sub Y{ return (-$_[1]+$r*$_[0]-$_[0]*$_[2]); } sub Z{ return (-$b*$_[2]+$_[0]*$_[1]); } # initial condition $x=1.0; $y=1.0; $z=1.0; print "$x $y $z \n"; # time step $h=0.01; $N=2000; # time development for($j=0; $j<$N; $j++){ $x1=$h*X($x,$y,$z); $y1=$h*Y($x,$y,$z); $z1=$h*Z($x,$y,$z); $x2=$h*X($x+$x1/2,$y+$y1/2,$z+$z1/2); $y2=$h*Y($x+$x1/2,$y+$y1/2,$z+$z1/2); $z2=$h*Z($x+$x1/2,$y+$y1/2,$z+$z1/2); $x3=$h*X($x+$x2/2,$y+$y2/2,$z+$z2/2); $y3=$h*Y($x+$x2/2,$y+$y2/2,$z+$z2/2); $z3=$h*Z($x+$x2/2,$y+$y2/2,$z+$z2/2); $x4=$h*X($x+$x3,$y+$y3,$z+$z3); $y4=$h*Y($x+$x3,$y+$y3,$z+$z3); $z4=$h*Z($x+$x3,$y+$y3,$z+$z3); $x+=($x1+2*$x2+2*$x3+$x4)/6; $y+=($y1+2*$y2+2*$y3+$y4)/6; $z+=($z1+2*$z2+2*$z3+$z4)/6; print "$x $y $z \n"; }

Ruby: # lorenz.rb def X(x, y, z) return (-Sigma*(x-y)) end def Y(x, y, z) return (-y+R*x-x*z) end def Z(x, y, z) return (-B*z+x*y) end Sigma, R, B=10.0, 28.0, 8.0/3.0 # 定数は大文字を使う！ # initial condition x, y, z=1.0, 1.0, 1.0 print(x," ",y," ",z,"\n") # time step h=0.01 N=2000 # time development for j in 1..N x1=h*X(x,y,z) y1=h*Y(x,y,z) z1=h*Z(x,y,z) x2=h*X(x+x1/2,y+y1/2,z+z1/2) y2=h*Y(x+x1/2,y+y1/2,z+z1/2) z2=h*Z(x+x1/2,y+y1/2,z+z1/2) x3=h*X(x+x2/2,y+y2/2,z+z2/2) y3=h*Y(x+x2/2,y+y2/2,z+z2/2) z3=h*Z(x+x2/2,y+y2/2,z+z2/2) x4=h*X(x+x3,y+y3,z+z3) y4=h*Y(x+x3,y+y3,z+z3) z4=h*Z(x+x3,y+y3,z+z3) x+=(x1+2*x2+2*x3+x4)/6 y+=(y1+2*y2+2*y3+y4)/6 z+=(z1+2*z2+2*z3+z4)/6 print(x," ",y," ",z,"\n") end Python: #lorenz.py

def X(x, y, z): return (-sigma*(x-y)) def Y(x, y, z): return (-y+r*x-x*z) def Z(x, y, z): return (-b*z+x*y) sigma=10.0 r=28.0 b=8.0/3.0 # initial condition x=1.0 y=1.0 z=1.0 print x,’ ’,y,’ ’,z # time step h=0.01 N=2000 # time development for j in range(N): x1=h*X(x,y,z) y1=h*Y(x,y,z) z1=h*Z(x,y,z) x2=h*X(x+x1/2,y+y1/2,z+z1/2) y2=h*Y(x+x1/2,y+y1/2,z+z1/2) z2=h*Z(x+x1/2,y+y1/2,z+z1/2) x3=h*X(x+x2/2,y+y2/2,z+z2/2) y3=h*Y(x+x2/2,y+y2/2,z+z2/2) z3=h*Z(x+x2/2,y+y2/2,z+z2/2) x4=h*X(x+x3,y+y3,z+z3) y4=h*Y(x+x3,y+y3,z+z3) z4=h*Z(x+x3,y+y3,z+z3) x+=(x1+2*x2+2*x3+x4)/6 y+=(y1+2*y2+2*y3+y4)/6 z+=(z1+2*z2+2*z3+z4)/6 print x,’ ’,y,’ ’,z [類題５] 非線形バネ 非線形バネの運動方程式 d2x dt2 =−ω 2 0x− βx3 を適当な初期条件のもとで数値的に解くには，新しい変数p = dx/dtを導入して，１階連

立の微分方程式に直して例題の方法を適用すればよい． dx dt = p, dp dt =−ω 2 0x− βx3≡ f(x) ここでは「リープ・フロッグ法」による解を示す．すなわち p_j+1 2 = pj+ h 2 · f(xj) xj+1 = xj+ h· pj+1 2 pj+1 = pj+1 2 + h 2 · f(xj+1) のように，間に p_j+1 2 をはさんで（蛙跳び），(xj, pj)→ (xj+1, pj+1) の１ステップを進め るのである．

6

行列計算

[例題６] ポアソン方程式 １次元ポアソン方程式 d2φ/dx2=−ρを差分化すると φj+1− 2φj+ φj−1 h2 =−ρj であるから，行列形式に直して         2 −1 0 · · · 0 −1 2 −1 · · · 0 · · · · · · · · 1 2 −1 0 · · · 0 −1 2                 φ1 φ2 · · · φN−1 φN         =         h2_ρ 1+ φ0 h2ρ2 · · · h2ρ_N−1 h2ρN + φN +1         となる．φ0, φN +1は境界条件として指定されるから，電荷分布ρj が与えられれば，右辺 のベクトルは既知となる．よって，この行列方程式を消去法で解けば，電位φj が求まる． 左辺のような対角要素とその両隣りだけがゼロでない行列を３重対角行列といい，物理の いろいろな場面に登場し，その数値解法もよく知られている． C: /* poisson.c */ #include <stdio.h> #define N 100

/* tri-diagonal matrix solver */

gauss(double diag[], double sub[], double sup[], double b[], int M){ int j; double t; for(j=0; j<(M-1); j++){ t=sub[j]/diag[j]; diag[j+1] -=t*sup[j]; b[j+1] -=t*b[j]; } b[M-1]/=diag[M-1]; for(j=(M-2); j>=0; j--){ b[j]=(b[j]-sup[j]*b[j+1])/diag[j]; } } main(){ int j; double L=1.0, h;

double diag[N], sup[N-1], sub[N-1], b[N]; double phi[N+2], rho[N+2], rho0=1.0;

/* boundary condition */ phi[0]=0; phi[N+1]=0; h=L/(double)(N+1); for(j=1; j<=N; j++){ rho[j]=rho0; b[j-1]=h*h*rho[j]; } /* N-th order matrix */ for(j=0; j<N; j++){ diag[j]=2; } for(j=0; j<(N-1); j++){ sub[j]=-1; sup[j]=-1; } /* solve */

gauss(diag, sub, sup, b, N); /* output */ for(j=1; j<=N; j++){ phi[j]=b[j-1]; } for(j=0; j<=(N+1); j++){ printf("%d %lf\n", j, phi[j]); } } Java: // poisson.java Perl: # poisson.pl Ruby: #poisson.rb Python: # poisson.py

from numpy import * def coef_matrix(n):

A=zeros([n, n], float) for j in range(n):

for j in range(n-1): A[j, j+1]=-1.0 A[j+1, j]=-1.0 return A # model 0, (1, 2, ... , N), N+1 N=100 L=1.0 # length h=L/(N+1) # initialization

phi=zeros(N+2, float) # potential rho=zeros(N, float) # charge density for j in range(N): rho[j]=1.0 # solve C=coef_matrix(N) K=linalg.inv(C) rhs=h*h*rho rhs[0]+=phi[0] rhs[N-1]+=phi[N+1] solution=dot(K, rhs) # output for j in range(1, N+1): phi[j]=solution[j-1] for j in range(N+2): x=h*j # position print x,’ ’,phi[j]