C 言語第 6 回 1 数値シミュレーション :2 階の微分方程式 ( シラバス10 11 回目 ) 1 2 階の微分方程式と差分方程式微分方程式を 2 d x dx + c = f ( x, t) 2 dt dt とするこれを 2 つの 1 階の微分方程式に変更する ìdx = y 2 2 d

Size: px

Start display at page:

Download "C 言語第 6 回 1 数値シミュレーション :2 階の微分方程式 ( シラバス10 11 回目 ) 1 2 階の微分方程式と差分方程式微分方程式を 2 d x dx + c = f ( x, t) 2 dt dt とするこれを 2 つの 1 階の微分方程式に変更する ìdx = y 2 2 d"

さやてらわ
5 years ago
Views:

1 C 言語第 6 回 1 数値シミュレーション : 階の微分方程式 ( シラバス10 11 回目 ) 1 階の微分方程式と差分方程式微分方程式を d x dx + c = f ( x, t) とするこれをつの 1 階の微分方程式に変更する ìdx = y d x dx d x dx ï dt c f ( x, t) c f ( x, t) + = Þ = - + Þ í ï dy = - cy + f ( x, t) = F( x, y, t) ïî dt それぞれの 1 階の微分方程式に Runge Kutta 法を適用するそれぞれ差分に直してより dx x - x dy y - y Þ, Þ dt Dt dt Dt x x n+ 1 n n+ 1 n k + k + k 6 + k k + k + k 6 + k n+ 1 - n = y を用いるここで k1 = f ( xn, yn, tn) Dt n+ 1 - yn = D k D t k t t k f x y + D t + D 1 1 = ( n +, n, n ) k D t k t t k f x y + D t + D k = f ( x + k D t, y + k D t, t + Dt) 3 = ( n +, n, n ) 4 n 3 n 3 n t RLC 回路のシミュレーション S E L R C x S E L R C x

2 C 言語第 6 回図のようなRC 回路にコイル (L) が入っているスイッチSで端子を電池 (E) 側から切り替えるとコンデンサに生じた電位が徐々に消滅していくそのときのコンデンサの電圧変化 x は時間とともに変化する電流の向きを RからLの方向にとると d x dx + + = 0 LC RC x となるこれをつの1 階の方程式にするには ì R c = d x dx d x R dx x ï L LC + RC + x = 0 Þ + = - Þ í dt L dt LC ï x f ( x, t) = - ïî LC なのでになる初期条件としてにする dx = y = f1 ( x, y) dt dy dy R 1 RCy + x = - cy + f ( x, t) Þ = - y - x = - = f x, y L LC LC x ( 0) ( ) = E y 0 = 0 ( ) 3 数値計算に必要な工夫微分方程式は ìdx = y ï dt í ï dy RCy + x = - ïî dt LC になるただし dy RCy + x = - はこのまま数値計算すると例えば dt LC のとき R=0.05(kΩ) C=0.5 (μf) L=10(mH) 3-6 dy RCy + x 0.05kW 0.5mF y + x W F y + x = - = - = dt LC 10mH 0.5mF H F y + x = = y + x = y + x となり係数が非常に大きな値となり数値計算には不適であるこのようなときには時間変数 t は秒で計測するのであるがこの 1 1 計測時間の単位を ms( = s ) や m s( = s ) のように短くする , 000

3 C 言語第 6 回 3 とよいこのことを見るために実際に数値を用いて dx 9 y= dt 4 10 dy æ ö d x æ dx ö = - ç y + x Þ = - ç x dt è.5 ø dt è dt.5 ø なので大きな数が現れるそこで 8 10 ( 後で使うため ( ) 4 10 を見越している ) で両辺を割ると d x æ 4 dx 10 ö d x æ dx 10 ö = ç + x Þ = - + x è ø 10 dt ç 10 dt 10.5 è ø æ 0.5 dx 10 ö = - ç + x 4 è 10 dt.5 ø になる分母の大きな数を時間 t に含めると æ ö æ ö ç x ç dt è dt ø 4 4 d ( 10 t) ç d ( 10 t) è ø d x 0.5 dx 10 d x dx 1 = - + x Þ = そこで新しい変数として 4 T = 10 t をとればもとの微分方程式は d x æ dx 1 ö = x T 10 t ç + = dt è dt 4 ø になり大きな数は姿を消す 4 ( ) 4 ことがわかるその結果プログラムでは T = 10 t について 0 秒から 60 秒間計測するとすれば t による実際の計測時間は 4 ( ) ì プログラムでの経過時間 : T = 60s T = 10 t T = 0sü ï t 0s ý Þ í = 1 60s 6 þ ï計測時間 : t = = s = î 6ms になる従って計測時間は ms で表示で表示するとよいこの変換を系統的に行うには dy RCy x d x R dx 1 = - Þ = - + dt LC dt L dt LC + æ ç x ö è ø -8 において LC = なので 10 8 のかわりに 1 LC を用いてと両辺を 1 LC で割る = 両辺に LC を掛ける

4 C 言語第 6 回 4 æ 1 ö æ ö = - ç + x Þ LC = - ç RC + x è ø è ø d x R dx d x dx dt L dt LC æ ö æ ö d x ç RC dx ç C dx Þ = - ç + x = - ç R + x æ t ö LC t L t d d d ç ç ç è LC ø è LC ø è LC ø そこで新しい変数として t T = LC をとればもとの微分方程式は d x æ C dx ö æ t ö 計測時間 t = LCT = - R x T : [ s] dt ç + L dt ç = プログラムでの経過時間 Þ è ø è LC ø msで表示 : LCT 10 3 になる従ってプログラムで用いる微分方程式は dx = y = f1 ( x, y) dt dy æ C ö = - R y + x = f ( x, y) dt ç L è ø になる 4 プログラム 1) ( ) ( ) (, ) ìd x = f1 x y Dt dx Dx ï = y = f1 x, y Þ = f1 x, y Þ í Dx Þ xnext = x + f1 ( x, y) Dt dt Dt ï îxnext - x = f1 ( x, y) Dt プログラムでは STEP æ ö xnext = x + f1 ( x, t) D t = calculate nextx x, y, Dt ç è ø ) R y x f ( x, y) f ( x, y) (, ) ìd y = f x y Dt dy æ C ö Dy ï = - Dy dt ç + L = Þ = Þ í Dt è ø ï î ynext - y = f ( x, y) Dt プログラムでは STEP æ ö ynext = y + f ( x, y) D t = calculate nexty x, y, Dt ç è ø ( ) Þ y = y + f x y Dt next, #include <stdio.h> #include <math.h> repo1.c #define R (0.05*1000.0) #define C (0.5/ ) #define L (10.0/1000.0)

5 C 言語第 6 回 5 #define INITIAL_TIME 0.0 #define LAST_TIME 10.0 #define INITIAL_X 1.0 #define INITIAL_Y 0.0 #define STEP 0.04 double f1(double x, double y) double fy; fy = y; return fy; double f(double x, double y) double fx; fx = -(R*sqrt(C/L)*y+x); return fx; double calclate_nextx(double x, double y, double step) double k1, k, k3, k4; double next_x; k1 = f1(x, y); k = f1(x+k1*step/.0, y+k1*step/.0); k3 = f1(x+k*step/.0, y+k*step/.0); k4 = f1(x+k3*step, y+k3*step); next_x = x + (k1+.0*k+.0*k3+k4) *step/6.0; return next_x; double calclate_nexty (double x, double y, double step) double k1, k, k3, k4; double next_y; k1 = f(x, y); k = f(x+k1*step/.0, y+k1*step/.0); k3 = f(x+k*step/.0, y+k*step/.0); k4 = f(x+k3*step, y+k3*step); next_y = y + (k1+.0*k+.0*k3+k4)*step/6.0; return next_y; void writedata(double t, double x) // ms で表示するため 1.0E3 を掛けておく printf("%7.4f,%10.3e\n", sqrt(l*c)*t*1.0e3, x); int main(void)

6 C 言語第 6 回 6 double t, last_t; double x, next_x, y, next_y; x = INITIAL_X; y = INITIAL_Y; t = INITIAL_TIME; last_t = LAST_TIME*(1.0+STEP); while(t < last_t) writedata(t, x); next_x = calclate_nextx(x, y, STEP); next_y = calclate_nexty(x, y, STEP); x = next_x; y = next_y; t = t+step; return 0; になる 5 表計算ソフトの自動立ち上げ(Windows 用のみ ) 計算結果を画面表示と同時にファイルを作成し作成したファイルを自動的に表計算ソフトで表にして表示するようにする 1)csv ファイルを作成 :fopen 関数 fopen( 番地, 決まった文字列 ); fopen(csv ファイル名の格納番地, "wt"); FILE * file_open; file_open= fopen(filename, "wt"); file_open の数値が零 (NULL) だと作成に失敗 )csv ファイルへの書き込み :fprintf 関数 fprintf(file * 型番地,...); printf("%5.f %10.3e\n", t, x); fprintf(file_open, "%5.f, %10.3e\n", t, x); 3) 書き込み終了後 csv ファイルを閉じる :fclose 関数 fclose(file * 型番地 ) fclose(file_open) 4) 表計算ソフトでcsv ファイルを開く :ShellExecute 関数 ShellExecute((HWND)0, "open", 番地, NULL, NULL, SW_SHOWNORMAL) ShellExecute((HWND)0, "open", csv ファイル名の格納番地, NULL, NULL, SW_SHOWNORMAL) 戻り値は特別で HWND 型になっている戻り値が3 以下だと csv ファイル用表計算ソフトの起動に失敗

7 C 言語第 6 回 7 の関数を使用する ShellExecute 関数を使用するために windows.h を #include することメニュー [ プロジェクト ]-[ プロパティー ] の [ 文字セット ] マルチバイト文字セットを使用するを選択プログラムは repo1.c に太文字部を追加や変更します #define _CRT_SECURE_NO_WARNINGS #include <windows.h> #include <stdio.h> #include <math.h> #define R (0.05*1000.0) #define C (0.5/ ) #define L (10.0/1000.0) #define INITIAL_TIME 0.0 #define LAST_TIME 50.0 #define INITIAL_X 1.0 #define INITIAL_Y 0.0 #define STEP double f1(double x, double y) double fy; fy = y; repo.c return fy; double f(double x, double y) double fx; fx = -(R*sqrt(C/L)*y+x); return fx; double calclate_nextx(double x, double y, double step) double k1, k, k3, k4; double next_x;

8 C 言語第 6 回 8 k1 = f1(x, y); k = f1(x+k1*step/.0, y+k1*step/.0); k3 = f1(x+k*step/.0, y+k*step/.0); k4 = f1(x+k3*step, y+k3*step); next_x = x + (k1+.0*k+.0*k3+k4)*step/6.0; return next_x; double calclate_nexty (double x, double y, double step) double k1, k, k3, k4; double next_y; k1 = f(x, y); k = f(x+k1*step/.0, y+k1*step/.0); k3 = f(x+k*step/.0, y+k*step/.0); k4 = f(x+k3*step, y+k3*step); next_y = y + (k1+.0*k+.0*k3+k4)*step/6.0; return next_y; void writedata(file *file_open, double t, double x) // ms で表示するため 1.0E3 を掛けておく printf("%7.4f,%10.3e\n", sqrt(l*c)*t*1.0e3, x); if(file_open!= NULL) fprintf(file_open, "%7.4f, %10.3e\n", sqrt(l*c)*t*1.0e3, x); int main(void) // 作成ファイル名は example.csv char filename[] = "example.csv"; // 作成ファイル追尾用変数 FILE * file_open = NULL; double t, last_t; double x, next_x, y, next_y; x = INITIAL_X; y = INITIAL_Y; t = INITIAL_TIME; last_t = LAST_TIME*(1.0+STEP); // ファイルの作成に失敗すると数値 (file_open=null) を返す file_open = fopen(filename, "wt"); writedata(file_open, t, x); while(t < last_t) next_x = calclate_nextx(x, y, STEP); next_y = calclate_nexty(x, y, STEP); x = next_x; y = next_y; t = t+step; writedata(file_open, t, x); // NULL のときにファイル作成失敗

$C 言語第 6 回 9 if(file_open == NULL) printf("\n====> ファイル (%s) の作成に失敗しています \n====> 既にファイルが開かれているかもしれません ", filename); getchar(); else HINSTANCE hinst; fclose(file_open); hinst =$

9 C 言語第 6 回 9 if(file_open == NULL) printf("\n====> ファイル (%s) の作成に失敗しています \n====> 既にファイルが開かれているかもしれません ", filename); getchar(); else HINSTANCE hinst; fclose(file_open); hinst = ShellExecute((HWND)0, "open", filename, NULL, NULL, SW_SHOWNORMAL); if(hinst <= (HINSTANCE)3) printf("\n====> 表計算ソフトの自動起動に失敗しました "); getchar(); return 0; エクセルからグラフ作成方法散布図 - 散布図 ( 平滑線 )

10 C 言語第 6 回 10 グラフの移動新しいシート

11 C 言語第 6 回 11 完了

C 言語第 6 回 1 第 6 回目レポートレポートの回数 (6 回目 ) 学生証番号と氏名を明記すること 1) RLC 回路においての点を提出必ず表紙を付け最後のページ ( の添付用 ) を表紙の次に入れる (A4 レポート用紙使用のこと ) d x dx コンデンサー電圧の満たす方程式 LC + RC + x = 0 を導けキルヒホッフの第 Ⅱ 法則を用いる

12 C 言語第 6 回 1 第 6 回目レポートレポートの回数 (6 回目 ) 学生証番号と氏名を明記すること 1) RLC 回路においての点を提出必ず表紙を付け最後のページ ( の添付用 ) を表紙の次に入れる (A4 レポート用紙使用のこと ) d x dx コンデンサー電圧の満たす方程式 LC + RC + x = 0 を導けキルヒホッフの第 Ⅱ 法則を用いるコンデンサーに蓄えられる電荷を Q とするとコンデンサの電圧 V ( x) 電流は Q を使って書ける = は Q を使って書ける Q の微分方程式をコンデンサの電圧 V ( x) = の微分方程式に変えて初期条件 y ( 0) = 0 は何を意味するか物理の言葉で説明せよ ) repo. のプログラムを作成し t-x グラフを作成せよプログラム t-x グラフ ( 散布図を選択 ) ができるが数値などの見栄えを美しくする事 ) グラフが示すように振動しながら減衰してゆくが方程式を元にして何故かの説明の3 点を提出 3) m = 0.1 kg の質点を先端にもつバネの減衰振動の方程式 d x dx m = -kx - b d x dx を LC + RC + x = 0 のプログラムを参考にして作成し数値 ìk = 0.1 A) í îb = 0.01

13 C 言語第 6 回 13 ìk = 0.1 B) í îb = 0.4 を用いて t-x グラフを作成せよの 3 点を提出プログラム :A) と B) 共通で 1 点時間間隔は秒 =5ms に #define STEP ただしデータの表示時間は t でよい : printf("%5.f,%10.3e\n", t, x); t-x グラフ :A) と B) の点 A) のグラフ B) のグラフ

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

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

16 C 言語第 6 回 16 第 6 回目レポート ( 添付用 ) d x dx 1) RLC 回路のコンデンサー電圧の満たす方程式 LC + RC + x = 0 を導け ) repo.1c/repo. のどちらかのプログラムを作成し t-x グラフを作成せよプログラム実行結果 t-x グラフグラフが示すように振動しながら減衰してゆくが方程式を元にして何故かの説明の4 点を提出 3) m = 0.1 kg の質点を先端にもつバネの減衰振動の方程式 d x dx m = -kx - b d x dx を LC + RC + x = 0 のプログラムを参考にして作成し数値 ìk = 0.05 A) í îb = 0.0 ìk = 0.1 B) í îb = 0.4 を用いて t-x グラフを作成せよ #define STEP 0.1 とし 50 秒間計測するの 5 点を提出プログラム :A) と B) 共通で 1 点実行結果 :A) と B) の点 t-x グラフ :A) と B) の点

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 ( C 言語第 8 回複素微分方程式の解法複素数の係数を持つ階の微分方程式複素数を z として微分方程式は dz dt = である特にとする f ( z, t) ( ) 実際にはが含まれていないので ( ) f ( z, t) = - 0.+ 0.5 i z Ü t f ( z) = - 0.+ 0.5i z Ruge-Kutta( ルンゲクッタ ) 法解くべき差分方程式はオイラー法で

C 言語第 6 回 1 数値シミュレーション :2 階の微分方程式 ( シラバス10 11 回目 ) 1 2 階の微分方程式と差分方程式微分方程式を 2 d x dx + c = f ( x, t) 2 dt dt とする これを 2 つの 1 階の微分方程式に変更する ìdx = y 2 2 d

C 言語第 6 回 1 数値シミュレーション :2 階の微分方程式 ( シラバス10 11 回目 ) 1 2 階の微分方程式と差分方程式微分方程式を 2 d x dx + c = f ( x, t) 2 dt dt とするこれを 2 つの 1 階の微分方程式に変更する ìdx = y 2 2 d