C 言語第 8 回複素微分方程式の解法 1 1 複素数の係数を持つ 1 階の微分方程式複素数を z として微分方程式は dz dt = である特にとする f ( z, t) ( ) 実際にはが含まれていないので ( ) f ( z, t) = i z Ü t f (

Size: px

Start display at page:

Download "C 言語第 8 回複素微分方程式の解法 1 1 複素数の係数を持つ 1 階の微分方程式複素数を z として微分方程式は dz dt = である特にとする f ( z, t) ( ) 実際にはが含まれていないので ( ) f ( z, t) = i z Ü t f ("

いおりままだ
5 years ago
Views:

1 C 言語第 8 回複素微分方程式の解法複素数の係数を持つ階の微分方程式複素数を z として微分方程式は dz dt = である特にとする f ( z, t) ( ) 実際にはが含まれていないので ( ) f ( z, t) = i z Ü t f ( z) = i z Ruge-Kutta( ルンゲクッタ ) 法解くべき差分方程式はオイラー法で k z - z = + f ( z, t) D t である第 5 回数値シミュレーション : 階の微分方程式 ( シラバス 8 9 回目 ) - 4 C 言語フロクラム :Ruge-Kutta 法では記号 k を用いて k + k + k + k 3 4 z+ - z = と変更されるここで k = f ( z, t) である k k f z t = ( +, k k f z t k = f ( z + k D t, t + ) 3 = ( +, 複素数操作プログラムの準備新しいプロジェクト bibu( 微分 ) を作成し以下のつの複素数を操作するプログラムをプロジェクトに追加する ) 以下のプログラムを [fukusosuu.h] として保存 : typedef struct tagcomplex double r; double i; COMPLEX; COMPLEX add(complex z, COMPLEX z); COMPLEX subtract(complex z, COMPLEX z); COMPLEX multiply(complex z, COMPLEX z); COMPLEX divide(complex z, COMPLEX z); COMPLEX cojugate(complex z); double absolute(complex z); COMPLEX multiply_umber(double x, COMPLEX z); Copyright Masaki Yasue, Dept. of Phys., Tokai Uiv, All rights reserved.

2 C 言語第 8 回 COMPLEX cexp(double re, double im); ) 以下のプログラムを [fukusosuu.c] として保存 : 該当箇所を第 7 回複素数の使用法 ( シラバス回目 ) のプログラムからコピー可能 #iclude <stdio.h> #iclude <math.h> #iclude "fukusosuu.h" // z+z COMPLEX add(complex z, COMPLEX z) COMPLEX z; z.r = z.r+z.r; z.i = z.i+z.i; retur z; // z-z COMPLEX subtract(complex z, COMPLEX z) COMPLEX z; z.r = z.r-z.r; z.i = z.i-z.i; retur z; // z*z COMPLEX multiply(complex z, COMPLEX z) COMPLEX z; z.r = z.r*z.r - z.i*z.i; z.i = z.r*z.i + z.i*z.r; retur z; // z/z COMPLEX divide(complex z, COMPLEX z) COMPLEX z; COMPLEX busi; if((z.r == 0) && (z.i == 0)) COMPLEX zero = 0, 0; pritf(" エラー :0 で割っています \"); retur zero;

3 C 言語第 8 回 3 busi.r = z.r*z.r + z.i*z.i; busi.i = z.r*z.i - z.i*z.r; z.r = busi.r/(z.r*z.r+z.i*z.i); z.i = busi.i/(z.r*z.r+z.i*z.i); retur z; // z^* COMPLEX cojugate(complex z) z.i = -z.i; retur z; // z double absolute(complex z) retur sqrt(z.r*z.r + z.i*z.i); // x*z COMPLEX multiply_umber(double x, COMPLEX z) z.r = x*z.r; z.i = x*z.i; retur z; 4 複素数操作関数の使用法複素数の操作を k + k + k + k 3 4 z+ - z = を例に取り解説するプログラム内ではに進む時時刻と計算して t 時刻 t ( t t ) + = + D k + k + k + k 3 4 z+ = z + 時刻 t ( t t) + = + D での z + を求めるここで k k = f ( z +, t = f ( z + k, t についてプログラムしてみよう : 引数は複素数と実数

4 C 言語第 8 回 4 f z 複素数実数 k D t D + + t Þ f(complex z, double t) (, t ) 計算結果は複素数 COMPLEX f(complex z, double t) k 複素数 z + = z + kの作り方 step step ì k double COMPLEX COMPLEX ü ï æ step ö w w multyply umber, k t k = 実数複素数 = = D ï ç z + k ï è ø ý Þ í w ï ï z z z, k k, D t stepþ ï k z + add ( z, w) îï k + k + k3 + k4 ルンゲクッタ法のプログラム部分 z + は関数 :calclate_extとして与えられる引数は変数 z 時間 t と時間の刻み calclate_ext(complex z, double t, double step) 計算結果は複素数 COMPLEX calclate_ext(complex z, double t, double step) k + k + k3 + k4 calclate_ext 関数 : z+ = z + ext z k + k + k3 + k4 z+ = z t + D を計算し計算結果は複素数 (COMPLEX) になり retur を使って計算結果を返す COMPLEX calclate_ext(complex z, double t, double step) COMPLEX k, k, k3, k4; COMPLEX ext_z; COMPLEX w, ww; k = f(z, t); Ü k = f ( z, t ) k // k = f ( z +, t = f ( z + k, t step/.0 step/.0 k k t w = multiply_umber(step/.0, k); Ü w= k : k = f ( z + D k, t step/.0 k w = add(z, w); Ü z + w : k = f ( z + k, t

5 C 言語第 8 回 5 w step/. 0 k k = f(w, t+step/.0); Ü k = f ( w, t : k = f ( z + k, t k // k3 = f ( z +, t = f ( z + k, t w = multiply_umber(step/.0, k); w = add(z, w); k3 = f(w, t+step/.0); // k = f ( z + k D t, t + D t) = f ( z + D tk, t + ) w = multiply_umber(step, k3); w = add(z, w); k4 = f(w, t+step); k + k + k3 + k4 w = multiply_umber(.0, k); Ü z+ - z = ww w k + k + k3 + k4 ww = add(k, w); Ü z+ - z = k + k + k3 + k4 w = multiply_umber(.0, k3); Ü z+ - z = ww = add(ww, w); Ü z ww w k + k + k3 + k - z = 4 + w ww k 4 k + k + k3 + k4 w = add(ww, k4); Ü z+ - z = w tep s w k + k + k3 + k4 + z D = k + k + k3 + 4 Ü z - = t k ww = multiply_umber(step/.0, w); ww ext z z D t ext_z = add(z, ww); Ü z+ = z + ( k + k + k3 + k4 ) retur ext_z; dz 微分方程式 f ( z, t) dt = は関数 :f(z,t): f ( z, t) ( ) = - + i z から与えられる

C 言語第 8 回 ( ) 引数は変数 z 時間 t f(complex z, double t) 計算結果は複素数 COMPLEX f(complex z, double t) f 関数 : f ( z, t) = (- 0.+ 0.

6 C 言語第 8 回 ( ) 引数は変数 z 時間 t f(complex z, double t) 計算結果は複素数 COMPLEX f(complex z, double t) f 関数 : f ( z, t) = ( ) f ( z, t) = i z を計算し計算結果は複素数 (COMPLEX) になり retur を使って返す COMPLEX f(complex z, double t) COMPLEX fz; a COMPLEX a = -0., 0.5; Ü i a Ü i z fz = multiply(a, z); ( ) retur fz; fz を省略した次の関数でもよい : COMPLEX f(complex z, double t) a COMPLEX a = -0., 0.5; Ü i a Ü i z retur multiply(a, z); ( ) 5 微分方程式解法プログラム z の初期値 ( t = 0 ) として COMPLEX z0 = 0, 0 ; ìt = 0 í Þ z 0 z = z0; î = t = 0; をとり 0. 秒おき ( D t = 0.) に t = 50秒まで計算する : D t = 0. Þ #defie STEP 0. t = 50 Þ #defie LAST TIME 50.0 計算結果をファイル名 :bibu.csv に書き込みエクセルを自動起動させるメニュー [ プロジェクト ]-[ プロパティー ] の [ 文字セット ] マルチバイト文字セットを使用するを選択 i z

7 C 言語第 8 回 7 以下のプログラムを [bibu.c] として先ほどのプロジェクト bibu( 微分 ) に追加する #defie _CRT_SECURE_NO_WARNINGS #iclude <widows.h> #iclude <stdio.h> #iclude <math.h> #iclude "fukusosuu.h" #defie LAST_TIME 50.0 #defie STEP COMPLEX f(complex z, double t) COMPLEX fz; COMPLEX a = -0., 0.5; fz = multiply(a, z); retur fz; COMPLEX calclate_ext(complex z, double t, double step) COMPLEX k, k, k3, k4; COMPLEX ext_z; COMPLEX w, ww; k = f(z, t); w = multiply_umber(step/.0, k); w = add(z, w); k = f(w, t); w = multiply_umber(step/.0, k); w = add(z, w); k3 = f(w, t); w = multiply_umber(step, k3); w = add(z, w); k4 = f(w, t); w = multiply_umber(.0, k); ww = add(k, w); w = multiply_umber(.0, k3); ww = add(ww, w); w = add(ww, k4); ww = multiply_umber(step/.0, w); ext_z = add(z, ww); retur ext_z;

8 C 言語第 8 回 8 void writetitle(file *file_ope) pritf("t, 実部,,t, 虚部 \"); if(file_ope!= NULL) fpritf(file_ope, "t, 実部,,t, 虚部 \"); void writedata(file *file_ope, COMPLEX z, double t) pritf("%5.f,%0.3e,,%5.f,%0.3e\", t, z.r, t, z.i); if(file_ope!= NULL) fpritf(file_ope, "%5.f,%0.3e,,%5.f,%0.3e\", t, z.r, t, z.i); it mai(void) // 作成ファイル名は bibu.csv char fileame[] = "bibu.csv"; // 作成ファイル追尾用変数 FILE * file_ope; double t, last_t; COMPLEX z, ext_z; COMPLEX z0 = 0, 0; z = z0; t = 0.0; last_t = LAST_TIME*(.0+STEP); // ファイルの作成に失敗すると数値 (file_ope=null) を返す file_ope = fope(fileame, "wt"); writetitle(file_ope); while(t < last_t) writedata(file_ope, z, t); ext_z = calclate_ext(z, t, STEP); z = ext_z; t = t+step; // NULL のときにファイル作成失敗 if(file_ope == NULL) pritf("\====> ファイル (%s) の作成に失敗しています \====> 既にファイルが開かれているかもしれません ", fileame); getchar(); else HINSTANCE hist; fclose(file_ope); hist = ShellExecute((HWND)0, "ope", fileame, NULL, NULL, SW_SHOWNORMAL); if(hist <= (HINSTANCE)3)

9 C 言語第 8 回 9 pritf("\====> 表計算ソフトの自動起動に失敗しました "); getchar(); retur 0;

10 C 言語第 8 回 0 第 8 回目レポートレポートの回数 (8 回目 ) 学生証番号と氏名を明記すること必ず表紙を付け最後のページ ( の添付用 ) を表紙の次に入れる (A4 レポート用紙使用のこと ) ) 3 つのプログラム (fukusosuu.h, fukusosuu.c, bibu.c) において 3 つプログラムプログラムを実行し計算結果データを出力して実部と虚部のつグラフ実部のグラフ虚部のグラフの 5 点ですグラフ作成時は縦軸の最小値を - 最大値を目盛間隔をにして表示する [ グラフの移動 ] を新しいシートにするを実行する事グラフは以下のようになります dz i z dt = - + を解き ) のグラフの振る舞いを説明せよ ) ( )

11 C 言語第 8 回 3) dz 0.05z 0.0iz dt = - + を解きプログラム : 自動作成のファイル名を repo.csv に変更すること // 作成ファイル名は repo.csv char fileame[] = "repo.csv"; ヒントは虚数 i は COMPLEX 型の 0, として定義する 0.0iz で z を作るには multiply 関数 iz を作るにはi と z に multiply 関数 0.0iz を作るには 0.0とiz に multiply_umber 関数を順に使う実部と虚部のつグラフ実部のグラフ虚部のグラフの 3 点を提出ですグラフ作成時は縦軸の最小値を 0 最大値を目盛間隔をにして表示する [ グラフの移動 ] を新しいシートにするを実行する事グラフは以下のようになります

12 C 言語第 8 回提出例 ( 使える表紙は次のページ ) コンヒュータ物理学演習 Ⅱ 第 8 回目月日次ページの添付用作成する解答など学籍番号氏名ページ目ページ目 3 ページ目以降

13 C 言語第 8 回 3 コンピュータ物理学演習 Ⅱ 第 8 回目月日提出学籍番号氏名

14 C 言語第 8 回 4 第 8 回目レポート ( 添付用 ) ) 3 つのプログラム (fukusosuu.h, fukusosuu.c, ryousi.c) において 3 つプログラムを作成プログラムを実行し計算結果 ( x y ( x, t) グラフ (_out = N_MAX/4 の箇所を変更します ) を作成 =0 =N_MAX/4 =N_MAX/ =3*N_MAX/4 =N_MAX の合計 7 点ですグラフ作成時は縦軸の最大値を 0 にして表示する [ グラフの移動 ] を新しいシートにするを実行する事 - ) データを出力して次のの場合の 5 つ ) グラフからは最初つの波束しかないがポテンシャルを通過するとつの波束に分かれてゆく様子が読み取れるつの量子をポテンシャルにぶつけるとつの量子に分裂し最後はつの量子が残るという事だろうか? この結果を説明せよ

C 言語第 7 回掛け算 (multiply number) ìz1 = x1 + iy1 í îz = x + iy 割り算 (devide number) ( )( ) ( ) Þ z z = x + iy x + iy = x x - y y + i y x + x y

C 言語第 7 回掛け算 (multiply number) ìz1 = x1 + iy1 í îz = x + iy 割り算 (devide number) ( )( ) ( ) Þ z z = x + iy x + iy = x x - y y + i y x + x y C 言語第 7 回複素数の使用法 ( シラバス 1 回目 ) 1 1 複素数複素数 (complex numbers) z は虚数単位 ìi í i = - î 1 を使ってつの実数 x, y から z = x + iy と作りますとくに x を z の実数部 (real part): x = Re( z) y を z の虚数部 (imarginary part): y = Im ( z)

C 言語第 8 回 複素微分方程式の解法 1 1 複素数の係数を持つ 1 階の微分方程式 複素数を z として 微分方程式は dz dt = である 特に とする f ( z, t) ( ) 実際には が含まれていないので ( ) f ( z, t) = i z Ü t f (

C 言語第 8 回複素微分方程式の解法 1 1 複素数の係数を持つ 1 階の微分方程式複素数を z として微分方程式は dz dt = である特にとする f ( z, t) ( ) 実際にはが含まれていないので ( ) f ( z, t) = i z Ü t f (