1 練習問題

1.1 内容

オイラー法で微分方程式を近似解を求めるプログラムを作成する．以下の練習問題を解け．

[練習1]

次の微分方程式をオイラー法で計算せよ．計算する$x$の範囲[0,2]とする．

$$\if 11 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi =2x x=0 y=0$$ (1)

• きざみ幅を0.1として，理論値と計算値を同じグラフに書け． 誤差のグラフも書け．
• きざみ幅を以下のように変化させた場合，$x=2$のときの誤差 を求めよ．

  0.1, 0.05, 0.01, 0.005, 0.001, 0.0005, 0.0001,

• きざみ幅と誤差の関係を考察せよ．
• きざみ幅が0.001の場合の誤差のグラフを描き，それを考 察せよ．横軸:$x$，縦軸:誤差．

[練習2]

同様にして，以下の微分方程式を計算せよ．

$$\if 11 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi =\cos x\sin y y(0)=0$$ (2)

[練習3]

以下の2階の微分方程式をオイラー法で計算せよ．

$$\if 12 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{2}... ...m{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi +y=\sin x y^\prime(0)=0 y(0)=0$$ (3)

1.2 ヒント

1.2.1 問1に関して

1.2.1.1 プログラム方法

アルゴリズムは，テキスト「常微分方程式の数値計算法」で説明しているので分かってい るものとする．まずは，計算の方針であるが，次のようにすればよろしい．

• 計算結果の$x_i$と$y_i$は，配列に格納する．すなわち，プログラム中では， $x_i$はx[i]に，$y_i$はy[i]という配列で表現する．
• 計算結果はグラフにするため，ファイルに保存する．第一列目:$x_i$の値，第二 列目:$y_i$の値，第三列目:誤差
• 誤差は，理論値から計算値$y_i$を引いたものとする．

この方針に基づいて，オイラー法の計算アルゴリズムを実現するフローチャートは図 1のようになる．いろいろな方法があるが，諸君にとって分かり やすく書いたつもりである．自分でプログラムが考えられる者は，これと同一にする必要 はまったく無い．

計算結果--$x_i$と$y_i$，誤差--は，ファイルには次のフォーマットで書けば良い．デー タの間はタブ(\t)で区切る．

   0.000000000000000	   0.000000000000000	   0.000000000000000
   0.100000000000000	   0.000000000000000	   0.010000000000000
   0.200000000000000	   0.020000000000000	   0.020000000000000
   0.300000000000000	   0.060000000000000	   0.030000000000000

                     長いので，この辺は省略

   1.700000000000000	   2.720000000000000	   0.170000000000000
   1.800000000000000	   3.060000000000000	   0.180000000000000
   1.900000000000000	   3.420000000000000	   0.190000000000001
   2.000000000000000	   3.800000000000000	   0.200000000000000

1.2.1.2 グラフの出力方法

計算したデータをグラフに出力する方法を述べる．詳細は，以前の gnuplotの講義ノートを見 れば分かるだろう．

先に示した結果を格納しているファイル名を「result.txt」として，話を進める．グラフ 描画アプリケーションはgnuplotを使う．

• gunplotを起動させるためには，「gnuplot」とターミナルに打ち込む．
• 第一列目をx軸，第二列目をy軸として，グラフを描画するためには，gnuplot起動 画面で以下のようにタイプする．

  
  

  	plot "result.txt" using 1:2
  

  
  線でつなぐ場合は，以下のようにする．

  
  

  	plot "result.txt" using 1:2 with line
  

  
  
• 第一列目の$x_i$と第二列目の$y_i$，そして理論解の$x^2$を一つのグラフに描く 場合，次のようにする．

  
  

  	plot "result.txt" using 1:2, x**2
  

  
  
• ディスプレイ上ではなくて，グラフをファイルにしたい場合がある．例えば， windowsでよく使われるemf(拡張メタファイル)ファイルとして，グラフを出力す る場合，次のようにする．
  1. 出力先をemfにする．

     
     

     	set terminal emf
     

     
     
  2. ファイル名を決める．

     
     

     	set output "hogehoge.emf"
     

     
     
  3. 後は，「plot」コマンドを使って，ディスプレイ上に描画するのと 同様にすれば，図のファイルができあがる．
• 適当なワープロで図を取り込めば，文書ができあがる．linuxでは， StarSuit Writerというワープロが使える．デスクトップ上の「StarSuit 7 文書 ドキュメント」をダブルクリックすると起動する．

図 1: オイラー方で常微分方程式を計算するフローチャート．

1.2.2 問2に関して

問1を参考にして，自分で考えろ．

1.2.3 問3に関して

これは，1階の連立微分方程式に式を変形する．

$$\if 11 \frac{\mathrm{d}y}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} y}{\mathrm{d}x^{1}}\fi =z$$ (4)

$$\if 11 \frac{\mathrm{d}z}{\mathrm{d}x} \else \frac{\mathrm{d}^{1} z}{\mathrm{d}x^{1}}\fi =-z-y-\sin x$$ (5)

これから，オイラー法の漸化式は

$$x_{i+1}=x_i+\Delta x$$ (6)

$$y_{i+1}=y_i+z\Delta x$$ (7)

$$z_{i+1}=z_i+(-z-y-\sin x)\Delta x$$ (8)

なる．後は，芋づる式に計算するだけ．
ホームページ: Yamamoto's laboratory
著者: 山本昌志