Untitled Document

物理を題材にしたC　その２

Cプログラムの勉強その１では、プログラム中では変数が実数なのか、整数なのか注意が必要であることを学びました。さて次は、IF文を練習します。IF文は変数がある値を取るかかどうかとういう判断を実行する文です。

つぎの例を見て下さい．

---------------------------------------------

#include <stdio.h>

main()

{

float x;

printf(" input a value : ");

scanf("%f",&x);

　　if(x==0.0)

printf(" your value is equal 0 %f \n",x);

　if(x>0.0)

printf(" your value is greater than 0 %f \n",x);

else

printf(" your value is smaller than 0 %f \n",x);

}

---------------------------------------------

最初のif文if(x==0.0)では実数型の変数xが 0.0（実数のゼロ）と等しいかどうか判断しています。Cでは等しいかどうかは==とイクオール記号を二つ重ねます。それで、xがゼロに等しい場合； your value is equal 0　と画面に印刷します。そうでない場合は,IF文の中ではなにもしません。このため、つぎのif文が実行されて　xが　ゼロより大きいかどうか判断します。その場合；

your value is greater than 0　と画面に印刷します。そうでない場合else以下が実行され； your value is smaller than 0　と画面に表示されます。この方法では明らかに問題が有ることがわかりますか？このプログラムをCコンパイラーでコンパイルして、実行して、xに0.0を代入して結果を見て下さい。何か変ですね？！！！どうしてなおしたらよいでしょうか？

実数型変数xへの入力がゼロのときのみ起こる、まれな場合の問題動作です。しかしプログラムではこれは「バグ」と呼ばれます。さてどう改良すればよいでしょうか？if文が２個有るので、問題が生じていることがわかりますか？ですから、if文を１個にしましょう。elseの代わりにelseif(x==0) printf();

を使ってみて下さい。

次に、微分を計算機に解かせてみよう。

時刻の関数　位置　x(t) の時間に関する微分は速度であること知っていると思います。これをdx/dtと書いておきます。これを計算するのによく使う方法はある時間間隔ｈの間で平均の速さを計算し、この時間間隔をゼロに持って行った極限が微分で言うと所の速度です。すなわち

　v(t)=dx/dt= Lim((x(t+h)-x(t))/h)

ところが計算機は上のような微分はできません、できるのは「数値微分もどき」で、これを差分と呼びます。上の式でLimを取ってvと結びます。

v(t)=(x(t+h)-x(t))/h

この式を次のように解きなおします。　x(t+h)=x(t)+h*v(t)

それで、この式を次のように読み直します。時刻tの位置x(t) と速度かける時間間隔h 　v(t)*hの和は,

時刻t+hの位置 x(t+h)を計算すると。計算機上では、h　が十分に小さければ（なにに比べて小さいのかは、考えてください）良い精度で微分を代行することができます。

その上、加速度は速度の時間微分ですね。これに対しても上と同じことを適用します。すなわち、時刻の関数　速度　v(t) の時間に関する微分は加速度でこれをdv/dtと書いておきます。これを計算するのによく使う方法はある時間間隔ｈの間で平均の速さを計算し、この時間間隔をゼロに持って行った極限が微分で言うと所の加速度a(t)です。すなわち

　a(t)=dv/dt= Lim((v(t+h)-v(t))/h)

差分方程式は

a(t)=(v(t+h)-v(t))/h

この式を次のように解きなおします。　v(t+h)=v(t)+h*a(t)

それで、この式を次のように読み直します。時刻tの速度v(t) と速度かける時間間隔h 　a(t)*hの和は,

時刻t+hの速度 v(t+h)を計算すると。

加速度を与えれば運動が判るってのが、古典力学＝ニュートン力学です。その例をやってみましょう。

今地上で物体を上に投げあげる場合を考えます。上方向をx軸に定義します。物体に働く力は-mgです。

運動方程式は　　F=m(dv/dt)=-mg 　　　両辺からmを消去して、

dv/dt=-g, 差分方程式では、v(t+h)=v(t)-g*h, x(t+h)=x(t)+v(t)*h

となります。これをCのプログラムでかくと（きれいではありませんが）

==========================================

/* solve differential equations **/
/* case of fall in the air */
#include<stdio.h>
main()
{
int n,nmax;
float t,x,v,g;
float dt,x1,v1;
nmax=1000 ;
g=9.8; /* gravitational const. */
dt=0.01 ; /* time step */
t=0.0,x=0.0,v=0.0;
/* initial values for the position and velocity */
for(n=0;n<nmax;n++) /* loop for the steps up to nmax-steps */
{
v1=v-g*dt; /* diffential equation of velocity in x */
x1=x+dt*v; /* differential equation of position in x */

t=t+dt;
printf(" %d %f %f %f \n",n,t,v,x);

v=v1; /* replace for the next step */
x=x1; /* replace for the next step */
}
}

==============================================

printf(" %d %f %f %f \n",n,t,v,x); 文で4つの量を毎回毎時間ごと印刷しているのは、gnuplotを使って、横軸時間t、縦軸速度vでプロットすれば、直線が得られるはずですね。また横軸時間t、縦軸位置xでプロットすれば二次関数がえられるはずです。これが見えればプログラムは正しく動いているという自信がつきます。　さて次にもうちょっと話をおもしろくします。空気抵抗を入れてみましょう。

空気抵抗力は　F(空気）= -k*m*v とかかれる場合を考えてみましょう。つまり速度に比例する空気抵抗力です。抵抗計数はkです。（質量をかけてプログラムが質量に無関係にしました）そうすると、

v1=v-g*dt-k*v*dt; /* diffential equation of velocity in x */

のようにたった一行変更すればよいですね！これで時間tと速度vをプロットすると、終端速度が見えますか？　　図の場合g=9.8m/s/s,k=0.1,dt=0.01で動かしてみました。10000ステップ！！

空気抵抗の入れ方は先にやった方法が唯一ではありません。もっと速度が速い場合、速度の２乗に比例することが知られています。これを今の例に無理矢理適用させてみましょう。ただし、問題は空気抵抗力が速度の１乗に比例する場合は、上向きの進行方向と下向きで速度の値は正から負に変化し運よく運動方程式を変更する必要がありませんでした。しかし、空気抵抗力が速度の２乗に比例する場合は、速度の正負は２乗されて恒に正になるため、外から人が方程式を変更させてやらねばなりません。これはif文でできますね。上向きに上がっているときは、

if(v>0.0)　 v1=v-(g+k*v*v)*dt;

とこんな風に書かねばなりません。下向きも変更してください。

　さて問題です。上向きの初速度をviとし、この物体が地上に落ちてきた時の速度をveとしましょう。また終端速度は計算もできますし、あるいは数値的に上の用に求めることもできます。この終端速度をvfとしましょう。有名な問題ですが、この３つの量には次の関係があります。

1/(ve*ve)=1/(vi*vi)+1/(vf*vf)

これを、数値的な差分を解く計算で示してみてください。できれば、地球以外の惑星でも、つまりgを変更して、やってみましょう。

===========================================================================================

次のプログラム上の技は、アレイと呼ばれるものです。これは１次元のアレイは物理ではベクトルと呼ばれますし、２次元のアレイは行列あるいはマトリックスと呼ばれます。３次元以上はテンソルですね。例えば１次元のアレイであるベクトルで4個の実数要素からかなるベクトル変数vは　float v[3] と記述されます。v[0],v[1],v[2] ,v[3]の各成分から成ります。同様に２次元の４個の各要素が実数の場合、　float w[3][3] のように書かれます。これはw[0][0],w[0][1],w[0][2],w[0][3],w[1][0],,,,,,,w[3][3]まで16個の成分から成ります。

ベクトルとマトリックスの積はベクトルですが、z=v・wのようには書けません。float z[3] に対して、　　z[0]=v[0]*w[0][0]+v[1]*w[0][1]+v[2]*w[0][2]+v[3]*w[0][3],

のように書き下す必要があります。これをfor文を使って書くこともできます。

for(k=0,k<4;k++){ for ( j=0;j<4;j++) { z[k] +=v[j]*w[k][j];　}}

皆さんはこれを使ってmatrix・vectorの計算をさせるプログラムを作ってください。

今回の課題として，コンデンサーの周りの静電ポテンシャルを計算することに挑戦してみましょう．ただし境界条件が決まっていないと問題は解けませんから，箱の周りは全てアースされていて電位つまりポテンシャル＝０，箱の中のコンデンサーの上側の板の電位を10，下側の電位を-10とします．このとき箱のなかの電位はどうなるだろうかというのが課題です．

電荷のないところでの電場Eはマックスウェルの方程式

div(E)=0

に従います．電場EとポテンシャルVの関係は

E=-grad(V)

ですから，

laplacian(V)=0

が解くべき方程式となります．これを具体的に書き下すと，

d2V(x,y)/dx2+d2V(x,y)/dy2=0

この方程式はラプラス方程式とよばれ，電磁気に限らずなぜか至る所にあらわれます．

さて計算機を使った数値計算では，微分を直接解くことはできませんから，差分で置き換えます．すなわち，xの関数f（x）の微分df/dxは，小さな量deltaを使って

{f(x+delta)-f(x)}/delta

と近似されます．

２次元の関数V(x,y)の場合も同様に，

dV(x,y)/dx={V(x+delta,y)-V(x,y)}/delta,

dV(x,y)/dx={V(x-delta,y)-V(x,y)}/delta,

２階の微分係数は１階の微分のもう１回の微分ですから，両者の差を取って，

d2V/dx2={dV(x,y)/dx-dV(x,y)/dx}/delta

={V(x+delta,y)-2V(x,y)+V(x-delta)}/delta2

y微分も同様に行い，結局ラプラス方程式の差分方程式は，

V(x+delta,y)-2V(x,y)+V(x-delta,y)+V(x,y+delta)-2V(x,y)+V(x,y-delta)=0

あるいはちょっと書き直して，

V(x,y)={V(x+delta,y)+V(x-delta,y)+V(x,y+delta)+V(x,y-delta)}/4

ラプラス方程式は，ある点の場の価Vは周りの上下左右４点の場の平均だと言っているのか！と感心してみてください．それでいろいろなところにラプラス方程式が顔をだすはずです．

さて，いよいよコンデンサーの周りの電位V(x,y)を求めてみましょう．適当な電位から出発して，ラプラス方程式の差分式をときながら，つまりV(x,y) を左右上下の価の平均値に置き換えて行きます．ただし，境界条件である，箱の部分とコンデンサーの電極の電位はそれぞれ一定に保っておきます．一度岳ではだめで，何度も繰り替えしを行うとだんだん望ましい調和のとれた配置の値に近付いて行くだろうと期待します．次にプログラムの例を示します．これもスタイルや変数の置き方の例であり，皆さんは皆さんのスタイルでプログラムを書いて下さい．

------------------------------------------

#include <stdio.h>

main()

{

double v[51][51],vl,vr,vu,vd;

int i,j,k,n,ix,iy,ly1,ly2,lx;

for( i=0;i<51;i++) /*　ポテンシャルの初期値と箱の電位*/

{ for (j=0;j<51;j++)

{ v[i][j] = 0.0;

}}

for (n=0;n<100;n++)　/* 繰り返す　*/

{ for(ix=1;ix<50;ix++) /*　箱の境界上は計算する必要なし*/

{ lx= (ix>=20 && ix<=30) ; /*コンデンサーのx条件 &&はAND論理積 */

for(iy=1;iy<50;iy++) /*　箱の境界上は計算する必要なし*/

{ ly1=(iy==25); /* コンデンサーのy条件*/

if(lx && ly1) /*論理積　AND*/

{

v[ix][iy]=10.0 ;　　/* コンデンサー上の電位*/

goto out1;　　　　/*　out1へ飛び出す*/

}

vl=v[ix-1][iy]; /* 左の点のポテンシャル*/

vr=v[ix+1][iy]; /* 右の点のポテンシャル*/

vu=v[ix ][iy+1]; /* 上の点のポテンシャル*/

vd=v[ix ][iy-1]; /* 下の点のポテンシャル*/

v[ix][iy]=(vl+vr+vu+vd)/4.0 ;　/* ラプラス方程式 */

out1:;

}

/* example print out */

ix=9 ;

for(iy=0;iy<51;iy++)

{

printf(" v(%d ,%d) = %lf \n",ix, iy,v[ix][iy]) ;

}

-------------------------------------------

このプログラム例では，ポテンシャルV(x,y)をひとつしか用意していないので，ラプラス方程式を利用して新しいV(x,y)を計算して格子点(x,y)をめぐってゆくとき，Vがどんどん置き換わって行きます．

最初のfor文で周囲の箱の電位を含めたすべての点の電位V(x,y)をv[i][j]=0.0として初期化していますが，皆さんはこの値をかえて任意の場Vから出発しても同じ結果にたどり着くことを確かめて下さい．またコンデンサーの電位を変化させて場がどのような形に変化するかも計算させてみて下さい．

３次元グラフを書かせることができたら、予想されるきれいな電位分布が見えるはずです。一つの方法としてはgnuplotを使ってみてはいかがでしょうか？x,y,V(x,Y)の３つの実数値を出力し、ファイルに保存し、これをgnuplotに渡して、splotしてみましょう。但しその前にset parametricとしてください。データをこちらが作った順番に読み込ませるためです。下にその例を示します。ここではコンデンサーの極板が１枚の例です。