[ ] π = C: /* pi.c $> gcc pi.c -lm $>./a.out / #include <stdio.h> #include <math.h> main(){ double pi=4atan(1); printf("%lf\n", pi); # pi=m

Size: px

Start display at page:

Download "[ ] π = 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"

ゆずさふじした
5 years ago
Views:

1 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! 1

2 [ ] π = 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: # pi.py $> python pi.py import math 2

3 print 4*math.atan(1) 3

4 2 [ ] F n+1 = F n + F n 1 F 0 = 1, F 1 = 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; 4

5 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 5

6 for j in 1..(n-1) ff=f0+f1 f0, f1=f1, ff return ff for k in print("f(",k,")=",fibo(k),"\n") 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! 6

7 [ ] 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 7

8 $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) return ff for j in print(j,"!=",factorial(j), "\n") 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) 8

9 3 [ ] f(x) = 0 x x 0 x n+1 = x n f(x n) f (x n ) n f(x ) = 0 x f(x) f(x) = x = 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); 9

10 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 a=2 x=a for n in x=rhs(x, a) p x 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 10

11 [ ] f(x) x 0 < x 1 f(x 0 ) f(x 1 ) < 0 x 0 < x < x 1 f(x) = 0 x f(x) = x 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; 11

12 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) a=2.0 left, right=0, a center=(left+right)/2; for n in

13 if residue(center, a)>0 then right=center else left=center center=(left+right)/2 p center 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 13

14 4 [ ] b a f(x)dx = N 1 j=0 h 2 (f(x j) + f(x j+1 )), ( h = b a ) N, x j = a + j h 0 x 1, f(x) = 1/(x 2 + 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++){ x=h*(double)j; 14

15 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)) 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) p s p Math.atan(1.0)) Python: # sekibun.py import math def f(x): 15

16 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 k 2 sin 2 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; 16

17 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{ 17

18 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)) N=100 pi=4.0*math.atan(1.0) h=(pi/2)/n for i in 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)) print(k, " ", s, "\n") 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 18

19 5 [ ] dx = σ(x y) dt σ = 10 dy = y + rx xz dt r = 28 dz = bz + xy dt b = 8/3 t = nh dx dt = f(x, t) k 1 = h f(x n, t), k 2 = h f(x n + k 1, t + h) x n+1 = x n (k 1 + k 2 ) k 1 = h f(x n, t), k 2 = h f(x n k 1, t + h 2 ), k 3 = h f(x n k 2, t + h 2 ), k 4 = h f(x n + k 3, t + h) x n+1 = x n (k k k 3 + k 4 ) 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)); double Y(double x, double y, double z){ 19

20 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: 20

21 // 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); 21

22 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"; 22

23 Ruby: # lorenz.rb 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, 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") Python: #lorenz.py 23

24 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 [ ] d 2 x dt 2 = ω2 0x βx 3 p = dx/dt 24

25 dx dt = p, dp dt = ω2 0x βx 3 f(x) p j+ 1 2 p j+ 1 2 = p j + h 2 f(x j) x j+1 = x j + h p j+ 1 2 p j+1 = p j+ 1 + h 2 2 f(x j+1) (x j, p j ) (x j+1, p j+1 ) 25

26 6 [ ] d 2 φ/dx 2 = ρ φ j+1 2φ j + φ j 1 h 2 = ρ j φ 1 φ 2 φ N 1 φ N = h 2 ρ 1 + φ 0 h 2 ρ 2 h 2 ρ N 1 h 2 ρ 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; 26

27 /* 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): A[j, j]=2.0 27

28 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] 28

Numerical Rosetta Stone 1 C, Java, Perl, Ruby, Python [ ] Hello world C: /* hello.c $> gcc hello.c $>./a.out */ #include <stdio.h> main(){ printf("hel

Numerical Rosetta Stone 1 C, Java, Perl, Ruby, Python [ ] Hello world C: /* hello.c $> gcc hello.c $>./a.out */ #include main(){ printf("hello world of C!\n"); Java: // hello.java $> javac hello.java