ドリフト運動（電場×磁場;垂直）（重力あり）

本日のソースコード

#include <>

main(){

double a,B,E,g,m,q,C,x,y,z,Vx,Vy,Vz,Vx1,Vy1,dt;

//C=q*B/m（定数）

int i,n;

a=0.0; //空気抵抗
B=2.0; //磁束密度//z方向
E=1.0; //電場//x方向
g=9.8; //重力加速度
m=1.0; //粒子質量
q=5.0; //粒子電荷
x=0.0; //
y=0.0; //
z=0.0; //
Vx=2.0; //
Vy=0.0; //
Vz=0.0; //


dt=0.0001; //時間（パラメータ）
for(i=0;i<100000;i++)
{
C=q*B/m;
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt+q*E*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
Vz=Vz-g*dt-a*Vz*dt;
printf("%e %e %e \n", x, y, z );
}
}

x-y平面への出力結果
x方向へ少しずつ広がっているのは微小時間tの値が大きいことによる誤差と考えられます。

x軸+方向に電場をかけた結果、ドリフト運動してy軸-方向へドリフト運動します。
さらに、振動の中心がすこしx軸+方向へずれているのが分かります。

z軸-方向へ重力場がはたらいているのでだんだん上昇速度が小さくなっています。

三次元plotその2 空気抵抗あり

右の図は下のプログラムの場合です


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

三次元のplot

先週のグラフを三次元で表してみました
（z軸方向に等速度運動）

gnuplot > plot "re"
→gnuplot > splot "re"
で三次元のplotが表示される
splotはファイルreに保存されたデータから一行目をx座標、二行目をy座標、三行目をz座標としてplotする。


#include

main(){

double a,x,y,z,Vx,Vy,Vz,Vx1,Vy1,C,dt;

int i,n;

a=0.0;
x=1.0;
y=0.0;
z=0.0;
Vx=2.0;
Vy=0.0;
Vz=0.1;
C=10;
dt=0.0001;
for(i=0;i<100000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;
z=z+Vz*dt;
Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e %e \n", x, y,z );
}
}

来週までの宿題

・プロットする点の減らし方
・磁場が変化する場合
・"三"次元プロット

磁場中の荷電粒子の運動（空気抵抗あり）

一様磁場：画面手前方向
電荷：＋
初速度：＋x軸方向

#include

main(){

double a,x,y,Vx,Vy,Vx1,Vy1,C,dt;

int i,n;

a=0.1;
x=1.0;
y=0.0;
Vx=2.0;
Vy=0.0;
C=10;
dt=0.0001;
for(i=0;i<1000000;i++)
{
x=x+Vx*dt;
y=y+Vy*dt;

Vx1=Vx+C*Vy*dt-a*Vx*dt;
Vy1=Vy-C*Vx*dt-a*Vy*dt;
Vx=Vx1;
Vy=Vy1;
printf("%e %e \n", x, y);
}

2009年10月26日月曜日へのコメント

return 0は関数の実行結果の値をゼロにするという意味だと思います
main関数でこうしてもあまり意味がないように思います

stdio.hは標準の入出力で利用する関数を呼び出すために必要な
情報を与えていると思えば良いと思います

来週の予告

磁場中の荷電粒子の運動（ローレンツ力あり）

運動方程式（xとyの連立方程式）
m*(dVx/dt)= q*Vy*B
m*(dVy/dt)=-q*Vx*B

どうやって解くか考えましょう。

重力場中（空気抵抗あり）の物体のx-tグラフのための数値計算

#include < stdio.h >
main()
{
double t, x, dt, dx, ta, tb, v;
double f(double ,double );
int i, n;
x = 0.0 ;
t = 0.0 ;
ta = 0.0 ;
tb = 100 ;
v = 0.0 ;
n = 1000;
dt = (tb - ta) / n ;
for ( i = 0; i < n ; i++){
x = x+ v*dt ;
v = v+ f(v, t)*dt ;
t = t + dt ;
printf(" %e %e \n", t, x);
}
}
double f( double x, double t){
double y, z, alpha, g;
alpha = 0.1;
g = - 9.8;
z = -alpha * x + g ;
return z;
}

来週のメニュー

来週は
2階の微分方程式（1次の項を含む）を
電解中の荷電粒子の運動（空気抵抗あり）をモデルに
やります。

ヒント
for()
{
x = x + v*dt
v = v + (A-B*v)*dt

}

(A=qE/m, B=alpha/m)

グラフを画像ファイルとして保存する

gnuplotで表示した画像の保存
clipboardに保存（グラフ画面の左上のアイコンをクリック）
GIMPを起動しファイルメニューから”画像の生成”を選び”クリップボードから”をクリック。するとグラフが表示される。名前をつけて保存を選び、形式を”JEPG画像”にして保存。

ブログにアップロードするには、ブログの投稿画面で”画像を追加”から”参照”

計算結果を保存する

コンパイルしたプログラムを
./a.out
で端末に計算結果が表示されるのでそれを
./a.out > ファイル名
で計算結果をファイルに保存する。

出来ました！

本当ですね。return 0がなくても実行できました。そもそもreturn 0とかstdio.hとかわけもわからずに使っているのですがこれらは一体何を意味しているんですか？



あとこのプログラムの実行は、最近、知り合いの人にダウンロードしてもらったcygwinでしています。残念ながら私にはダウンロードに関して、詳しい説明は無理ですがほかの二人なら自力でできるでしょう。是非ダウンロードしてみんなでどんどん進みましょう！！！

うまくいっているようですね

羽部です

うまくいっているようですね
ただ
return 0;
はいらないのでは

掛け算の表つくってみました

#include
int main(){
int i,j;
for(i=1;i<=9;i++){
for(j=1;j<=9;j++){
printf("%dx%d=%d\n",i,j,i*j);
} printf("****************\n");
}
return 0;←（いらない）

}

関数

プログラムの一部を別の関数にできる

そのためには

関数名をmain(){}で示して（宣言して）おく

関数のプログラムをつくる

データの渡し方は、
変数、ポインター、グローバル変数を使う、３つの方法がある

下に示したプログラムの例では変数でデータを関数に渡している
各自確認しよう

Eular法 のプログラム

#include < stdio.h >
main()
{
double t, x, dt, dx, ta, tb;
double f(double ,double );
int i, n;
x = 0.0 ;
t = 0.0 ;
ta = 0.0 ;
tb = 100 ;
n = 1000;
dt = (tb - ta) / n ;
for ( i = 0; i < n ; i++){
dx =f(x, t) *dt ;
x = x+ dx ;
t = t + dt ;
printf(" %e %e \n", t, x);
}
}
double f( double x, double t){
double y, z, alpha, g;
alpha = 0.1;
g = - 9.8;
z = -alpha * x + g ;
return z;
}

微分方程式の解き方

c言語の解説ホームページ

http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/

投稿の方法

右上の「投稿」をクリック
タイトルと本文を入力後、
下の「投稿を公開」をクリック

※注　他人も編集可能。

数値微分

前進差分で数値微分し，dxを半分ずつ小さくしたときの
数値微分値の変化を調べよ

微分する関数は
f(x)=x^4

同様のことを中心差分で行なえ

どれだけ大きい数の計算を行えるかを調べるプログラム

その１

#include < stdio.h >
main(){
int i=1,x;
for(x=1; x<70 ;x++){
i=i*2;
printf("%d %d \n" ,i,x);
}
}

その２

#include < stdio.h >
main(){
int count=1;
int x=1;
int y=0;
while (y !=x){
y=x;
x=x*2;
printf("%d\n",count,x);
count++;
}
}


while :繰り返しの記号
y !=x :y=xでないならば count :番号を付けていく xとyが一致しない限り繰り返されるプログラム。計算機が計算できる範囲を超えたらx=yとなるのでここで計算が止まる。

n項まで計算できるsinのテイラー展開のプログラムです

#include < stdio.h >
#include < math.h >
main(){
char FileName[124];
FILE *fp;

int i, j, k, n, m, e;
float a, b, c, d;
double x;
float dx, s;
x = 0.0e0;
dx = 0.1e0;
x = -dx;
b = 0;

printf("結果の保存先ファイル名を入力してください>");
scanf("%s" , FileName);
fp = fopen(FileName , "w");

printf("テイラー展開のn項まで計算します。nを入力してください>");
scanf("%d", &n);

for(e=1; e<=70; e++){
x = x + dx;

for(m=1; m<=n*2; m = m + 4){
b=0;
for(i=1; i<=m; i = i + 4){
a = 1.0;
for(j=1; j<=i; j++){
a = a * 1/j;
}
a = a * pow(x,i);
b = b + a;
}
}

for(m=3; m<=n*2; m = m + 4){
c = 0;
for(i=3; i<=m; i = i + 4){
a = 1.0;
for(j=1; j<=i; j++){
a = a * 1/j;
}
a = a * pow(x,i);
c = c + a;
}
c = -c;
}

d = b + c;
printf("%d %f %e\n", e, x, d);

fprintf(fp,"%d %f %e\n", e, x, d);
}

fclose(fp);

return 0;

}

実行したら、最初に実行結果を出力するファイル名を記入してください。
その後に、n項まで計算するための"n"を入力してください。

夏休み

８月３日が夏休み前の最後の日です
渡辺君も参加してください

cで扱える数のテスト

コンパイルの仕方

cのプログラムのコンパイルの仕方
コマンドは
gcc -lm
ここで -lmはmath.hを使うときにコンパイルのオプションとして
つける -lはリンク先を指定している
-lmは数学関数などのコンパイル済みのライブラリーを指定している

gcc test.c -lm

3-5 数値微分

微分値df/dxを数値的に求める(７月２７日の課題）
基本はテーラー展開


微分値の計算を右辺第一項で近似すると，


誤差はdxの一次に比例，これを一次精度と呼ぶ
これのプログラムを作ってみる
f(x)としては指数関数を使う

２次精度にするにはどうしたらよいか

3-4-3. gnuplotの使い方

使用例
ファイルreの２行目をx３行目をyとしてグラフを書く

astro-PBG4-12:~/tmp habe$ gnuplot <----- gnuplotの起動

G N U P L O T
Version 4.0 patchlevel 0
last modified Thu Apr 15 14:44:22 CEST 2004
System: Darwin 8.11.0

Terminal type set to 'x11'
gnuplot> plot "re"
gnuplot> plot "re" u 2:3
gnuplot> plot "re" u 2:3 with line
gnuplot> set yrange[-1:1]
gnuplot> plot "re" u 2:3 with line
gnuplot> exit

3-4-0. sin xのテーラー展開

3-4-2. プログラム例 （pow関数を使う）

#include < stdio.h >
#include < math.h>
main()
{
int i, n;
double x;
float dx, s ,c1,c3,c5,c7,c9;
x=0.0e0;
n= 70;
dx = 0.1e0;
x = -dx ;
c1=1;
c3=c1/2/3;
c5=c3 / 4 /5;
c7 = c5 /6 /7;
c9 = c7 /8 /9;
for (i=1;i<=n; i++ ){
x = x + dx;
s = x -c3 *pow(x,3) + c5*pow( x,5) - c7 * pow(x,7)+c9* pow(x,9);
printf(" %d %f %e \n", i , x, s);
}
}

3-4-1. プログラム例

#include < stdio.h >
main()
{
int i, n;
float x, dx, s ;
x=0.0;
n= 70;
dx = 0.1;
x = -dx ;
for (i=1;i<=n; i++ ){
x = x + dx;
s = x -1/6.0 *x*x*x +1/5.0/4.0/3.0/2.0* x*x*x*x*x
- 1/7.0/6.0/5.0/4.0/3.0/2.0* x*x*x*x*x*x*x ;
printf("%d %f %e \n", i , x, s);
}
}

3-4. テーラー展開で値を求めてみる

sin x のテーラー展開
n項までの和を求めるプログラム
結果をグラフにする
グラフを印刷してみる

3-3. 計算までの道筋

プログラムを書く
翻訳する（コンパイル）
実行する

3-2. c プログラムの例

これは何をしているのか,考えてみよう


#include
main()
{
int a ;
a = 1 + 2 ;
printf("%d¥n" , a);
a = 1 - 2 ;
printf("%d¥n" , a);
a = 1 * 2 ;
printf("%d¥n" , a);
a = 1.0 / 2 ;
printf("%d¥n" , a);
}

3-1. unix (Linux)の使い方

1) loginとlogout
他の計算機と接続するには sshを使う

2) ファイル

3) ディレクトリ
　ファイルをしまっておく場所で階層的に構成されている
　Linuxでは
　　一番大きな階層がルートディレクトリ(/)
　　その下のディレクトリの例(/bin, /usr, /homeなど）

4) 標準入出力とパイプラインとリダイレクト
　コマンドの出力は普通標準出力先（ディスプレイ）へ
　これをリダイレクトして例えばファイルに入力できる
　　ls > log.txt
これはlsの出力先を>でlog.txtというファイルに変更している

　コマンドの出力を別のコマンドの入力に出来る
　　ls > less　
これはlsの出力をlessに渡している
5) ファイル管理と複写、移動、削除，ディレクとの作成と削除，文字列の検索
ls, cp, mv, rm, mkdir, rmdir, grep

6) フォアグランドとバックグランド

3. unix (linux)

計算機の制御システムのこと
　c言語で書かれており、世界共同で開発され続けれている
計算を目的としている計算機で多く使われている
(pcからスーパーコンピューターまで）
複数の計算を同時に行なえる
複数の人が一つの計算機を同時に使用できる
解析ツールもそろっている(gcc, gnuplotなど）

2. プログラムとは

問題を計算機に計算させる手順を書いたもの，
ソースあるいはソースコードとも言う

使用する言語はc、c++あるいはfortran

取り上げるテーマ
惑星円盤における粒子の衝突による成長
膨張宇宙における密度揺らぎの成長

1. 内容（目次）

プログラムとは
Unix, Linuxについて
テーラー展開
運動方程式の解法

参加者

渡辺隼輔
海木寛之
佐藤祥太
若山真梨子
石倉未奈

羽部朝男

若山さん

若山さん
どうも、書き込めると思いますので
やってみてください

よろしく

羽部です

みなさん
連絡はメイルでしますが，
図や記録はここに投稿しておくと便利かと思い
作ってみました

また、セミナーで使ったファイルをのせていきますので
コメントください
よりよいものにしていきます

時間ができたら，整備します

関係者以外には公開されない
ブログになっていますので自由に書き込んでください