| PDFﾌｧｲﾙの印刷 |
|
【学内】PDFファイルリーダーで印刷可能です。PDFファイルでの配布のｐｄｆのリンクをクリックしてください。現われる講義録を で印刷してください。 |
 |
PDFファイル閲覧・印刷にはPDFファイルリーダーが必要です。 |
 |
無償で使用できる高性能なPDF作成ソフト |
|
|
|
|
|
|
第７回：数値シミュレーション（PDFファイルでの配布：c7.pdf、
微分法方程式の解法（補足）：csample.pdf）（Web版は更新されていません、pdfをご覧下さい） |
|
|
（教科書１８６ページの１３・５「微分方程式」参照）
目次
１．微分方程式と物理現象
２．微分方程式を差分方程式へ
３．Ｃ言語ﾌﾟﾛｸﾞﾗﾑ：Euler法
プログラム（repoS.c）
４．Runge-Kutta法
第七回目レポート
力学では、マニュアル書に従って「誰がやっても、過去の透視・未来の予言」ができます。時間を逆に辿ったり、先に進んだりします。これを使って透視や予言するのですが、道具は「水晶玉」等ではなくて、「時間についての微分方程式」になります。人間の能力では、実際の現象をあまりうまく予言できないので、スーパーコンピュータ等に予言してもらうことがしばしばです。台風の進路などの予言や、地球にぶつかりそうになる彗星の進路予言等、これですね。もっと身近には、１２：００発の新幹線に東京で乗ると大阪に何時につくか、予言できちゃってますね。
さて、ここでは微分方程式をコンピュータにといてもらうことをやってみます。使用する方法は
と、それを改良した
です。
微分方程式は、
である。 は微小変化分という意味で、微分では
とする。そこで、をやめて、適当な幅（微小幅）として近似する。また、
とかけるので、 がわかれば が計算できる。すなわち、最初の を与えて、次々に計算を進める。最初の値
を決めて を計算する。 を一定の値（0.01とか0.05とか0.1等）として、それを
とし、つぎの値を とすれば、
 、
となる。以下、同様にして
となる。そして、（ ）と並べていけば、答えが得られる。これを
という。ここで例として
dx/dt=x⇒x=et(=exp(t))
を解かせてみる。数値計算と答え：exp(t)を比較すれば・・・
このように、xが大きくなると誤差も少しずつ増えるが、大雑把にはあっている。
dx/dt=xをプラグラムする。
において
である。
| |
#include <stdio.h> |
| |
#include <math.h> |
| |
|
| |
#define STEP 0.05 |
| |
|
| |
// dx/dt=f(x)=x |
| |
double f(double x, double t) |
| |
{ |
| |
|
double fx; |
| |
|
|
| |
|
fx = x; |
| |
|
|
| |
|
return fx; |
| |
} |
| |
|
| |
// next_x - x = f(x,t) *step |
|
double
calclate_nextstep(double x, double t) |
| |
{ |
| |
|
double next_x; |
|
|
|
|
|
next_x = x + f(x,
t)*STEP; |
| |
|
|
|
|
return next_x; |
| |
} |
|
|
| |
int main(void) |
|
{ |
| |
|
double x, x_next, t; |
|
|
|
| |
|
x = 1.0; |
|
|
t = 0.0; |
| |
|
|
|
|
printf("%s, %s,
%s\n", "t", "x=f(t)", "exp(t)"); |
| |
|
|
|
|
while(t<(1.0+STEP)){ |
| |
|
|
printf("%f, %e, %e\n", t, x, exp(t)); |
| |
|
|
x_next = calclate_nextstep(x, t); |
| |
|
|
t = t+STEP; |
| |
|
|
x = x_next; |
| |
|
} |
| |
|
|
| |
|
getchar(); |
| |
|
|
| |
|
return 0; |
| |
} |
になる。
Euler法の改良版であり
の代わりに
を用いる。ここで、
である。これを図示すれば、すばらしい結果となる。
学生証番号と氏名を明記すること
必ず表紙を付け、このレポート問題を印刷し表紙の次にいれること
（Ａ４レポート用紙）
- １）repoS.cの
- A)プログラム
- B)main()プログラム内でtは０〜１の値を取らせている。そのためには「while(t<=1.0)」でよさそうであるが、「while(t<(1.0+step))」としてあるのは理由がある。その理由を述べよ。
- C)実行結果
- D)以下の画面表示をするように変更したプログラム：空白を使って調節しないこと。printfの書式設定を使用し「表題は左寄せ、数値は右寄せ」の設定にすること。
の３点を提出。
の２点を提出。
- ３）以下のRC回路でスイッチを入れた後のコンデンサCの電圧変化をRunge−Kutta法を用いて解くこと。但し、C=１μH、R=１ｋΩ、E=１Vとする。また、スイッチを入れた後の経過時間（t）をｍｓ（＝１０−３秒）で計るとする。コンデンサCの電圧はｘ＝f(t)とする。
の３点を提出。
|
