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 1 3 4 n+ 1 - n = y を用いる ここで k1 = f ( xn, yn, tn) Dt 1 3 4 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
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 0.05 10 W 0.5 10 F y + x = - = - = - -3-6 dt LC 10mH 0.5mF 10 10 H 0.5 10 F -3 9 9 0.05 0.5 10 y + x 6 10 4 10 = = 0.005 10 y + x = 0.5 10 y + x -9.5 10.5.5 となり 係数が非常に大きな値となり 数値計算には不適である このようなときには 時間変数 t は 秒で計測するのであるが この 1 1 計測時間の単位を ms( = s ) や m s( = s ) のように短くする 1000 1000, 000
C 言語第 6 回 3 とよい このことを見るために 実際に数値を用いて dx 9 y= 9 4 10 dt 4 10 dy æ ö d x æ dx ö = - ç 0.5 10 y + x Þ = - ç 0.5 10 + x dt è.5 ø dt è dt.5 ø なので 大きな数が現れる そこで 8 10 ( 後で使うため ( ) 4 10 を見越している ) で両辺を割る と 9 4 9 d x æ 4 dx 10 ö d x æ 0.5 10 dx 10 ö = - 0.5 10 ç + x Þ = - + x.5 8 8 8 è ø 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 Þ = - 0.5 + 8 4 10 10.5 4 そこで 新しい変数として 4 T = 10 t をとれば もとの微分方程式は d x æ dx 1 ö = - - 0.5 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 = = 10000 1000 s = î 6ms になる 従って 計測時間は ms で表示で表示するとよい この変換を系統的に行うには dy RCy x d x R dx 1 = - Þ = - + dt LC dt L dt LC + æ ç x ö è ø -8 において LC = 0.3 10 なので 10 8 のかわりに 1 LC を用いて と 両辺を 1 LC で割る = 両辺に LC を掛ける
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/1000000.0) #define L (10.0/1000.0)
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)
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 ファイル用表計算ソフトの起動に失敗
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/1000000.0) #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 0.005 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;
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 = ShellExecute((HWND)0, "open", filename, NULL, NULL, SW_SHOWNORMAL); if(hinst <= (HINSTANCE)3) printf("\n====> 表計算ソフトの自動起動に失敗しました "); getchar(); return 0; エクセルからグラフ作成方法 散布図 - 散布図 ( 平滑線 )
C 言語第 6 回 10 グラフの移動 新しいシート
C 言語第 6 回 11 完了
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
C 言語第 6 回 13 ìk = 0.1 B) í îb = 0.4 を用いて t-x グラフを作成せよ の 3 点を提出 プログラム :A) と B) 共通で 1 点 時間間隔は 0.005 秒 =5ms に #define STEP 0.005 ただし データの表示時間は t でよい : printf("%5.f,%10.3e\n", t, x); t-x グラフ :A) と B) の 点 A) のグラフ B) のグラフ
C 言語第 6 回 14 提出例 ( 使える表紙は次のページ ) コンヒ ュータ物理学演習 Ⅱ 第 6 回目 月 日 次ページの添付用 作成する解答など 学籍番号 氏名 1 ページ目 ページ目 3 ページ目以降
C 言語第 6 回 15 コンピュータ物理学演習 Ⅱ 第 6 回目 月日提出 学籍番号 氏名
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) の 点