. 演習課題 演習課題 ( 必修 ) 台形公式とシンプソンの公式を用いて次の定積分の値を求めよ S = + ( a 4 + b 3 + c + d + e) d なお, +, a, b, c, d, e は各自で設定せよ 教科書第 章を参照のこと y S は右図の + の面積に等しい この区間を 個の Δ に分割し を順に,,..., + それに対応する y を y, y,..., y + とおくと 面積 S は以下の公式により求めることができる y = f () 台形公式 0 + Δ Δ S = {( y + y ) + ( y + y3 ) + + ( y + y+ )} シンプソンの公式 Δ S = {( y + 4y + y3 ) + ( y3 + 4y4 + y5 ) + + ( y - + 4 y + y + )} 3 シンプソンの公式を使用する場合は 分割数は偶数であることに注意する [ 教科書のプログラムの変更 ] () 基本的には被積分関数が異なるだけで 教科書 章のプログラムを使用できる () これまでに扱っていない内容は Fucto プロシージャ ( 教科書 5-3) の使用である Fucto プロシージャを使う利点は メイン プログラムが汎用性を持ち Fucto プロシージャを変更するだけであらゆる関数に関する計算が可能になることである (3) できれば 教科書にある二つのプログラムを一つにまとめよう (4) a~e は Fucto プロシージャで直接指定するか またはセルから読み込むようにする セルから読み込む場合は a~e は次ページのように Fucto プロシージャに渡す必要がある (5) プログラムが文法上間違っていないのに結果がおかしい場合には 教科書にあるデバックの方法 (p.55) を利用してプログラムの修正を行う () プログラムリスト フローチャート 出力結果 公式や分割数の違いに関する考察をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 3
いくつかの引数を持つ Fucto プロシージャの形式 Fucto マクロ名 ( 引数 As データ型, 引数 As データ型, ) As データ型マクロ名 = 処理内容 ( 引数を使った計算式 ) Ed Fucto [ 使用例 ] BMI( ボディ マス インデックス ) の計算 解説 BMIcalc が Fucto プロシージャのマクロ名 IputBo で入力した身長と体重を (Double 型変数である ) Heght や Weght に格納した後 それを引数として BMI を計算する Fucto プロシージャに渡す Fucto プロシージャでは BMI 値の計算結果を (Double 型変数である )BMIcalc に格納し 呼び出し元の Sub プロシージャに戻り値として返す Sub プロシージャでは 戻り値を ( 変数 )BMI に格納し その内容を MsgBo で表示する Sub Prog() Heght# = IputBo(" 身長 (cm) を入力してください ") Heght# = Heght*0.0 Weght# = IputBo(" 体重 (kg) を入力してください ") BMI# = BMIcalc(Heght, Weght) If BMI < 8.5 The MsgBo "BMI は " & BMI & " やせすぎです " ElseIf BMI < 5 The MsgBo "BMI は " & BMI & " 適正体重です " ElseIf BMI < 30 The MsgBo "BMI は " & BMI & " 肥満です " Else MsgBo "BMI は " & BMI & " 高度肥満です " Ed If Ed sub Fucto BMIcalc ( Heght As Double, Weght As Double ) As Double '************************************************************************** ' 身長 (Heght(cm)) 体重(Weght(kg)) から BMI を計算する ' 引数 : Heght(Double), Weght(Double) ' 戻り値 : BMIcalc(Double) '************************************************************************** BMIcalc = Weght / (Heght*Heght) Ed Fucto 4
演習課題 ( 必修 ) 非線形方程式の二分法による解法 f()=e - という式で,f()=0 となる を求めよ ただし a = -0 b = 0 解の収束判定値(eps) = 0.0000 とする 教科書第 章を参照のこと [ 教科書のプログラムの変更 ] これまでに扱っていない内容は 複数の Sub プロシージャの使用 ( 教科書 5 章 5-) Et Sub(p.53) や Et For(p.40) の使用 である 課題との違い教科書では反復回数 (ma) を終えた時点の値を解としているが 課題では 解の候補 ( と 下図では a と b) の差の絶対値が十分小さく (0.0000 以下 ) なったときに解とする は大きめの値 ( 例えば 00) にしておこう または を更新した後に収束判定のための f 文を入れ 条件を満たしたときに解を出力し solve プロシージャから抜ければ (Et Sub) 良い Net の後の MsgBo では ma 回繰り返しても条件を満足しないことを出力しよう 注 ) p.97 Net の後の MsgBo にある Chr(3) & Chr(0) は vbcrlf( 課題集 p.7) と同じで改行を行うもの f() f() 新しい a f(b) f(a) a c f(c) c b f(a) f(b)<0 より a~b の間に f()=0 の点がある. c=(a+b)/ を計算し, f(c) を求めると f(a) f(c)>0, f(b) f(c)<0 より,c~b の間に解がある. そこで,c を a にして,[ 新しい a~b] 区間で同様に繰り返す. () プログラムリスト フローチャート 出力結果をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 5
演習課題 3 Naper の常数 e は 次式において として得ることができる e = +.7888845904535360 ( ) k = k! () k を入力データで与えて a k = = を計算せよ k! 3 ( k ) k () 上の結果を利用して e の値を求めよ (<50) は入力データとして与えよ フローチャートは以下のようになる ( スペースの関係で () は A とおいている ) () () 初め 初め k = 入力データ = 入力データ A b = s = 0 j = ~k の反復 = ~ の反復 b = b * j j の反復 A s = s + a 但し () の k は にする a = /b の反復 a を出力 e = + s 終わり e を出力 終わり () プログラムリスト 入力データ (k, の値 ) 出力結果 (a k, e の値 ) の値と e の精度の関係に関する考察をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 6
演習課題 4 地上から角度 a, 初速度 v 0 で発射された質点は 放物線を描いて再び地上に落下する この運動の軌跡を y 平面上に描け なお 角度 a, 初速度 v 0 は任意に設定してよいが 最大到達高度 移動距離等を計算して考察に加え プログラムの確認を行うこと 質点が発射されてから地上に到達するまでの時間 T をあらかじめ求めておく 時刻 t における質点の 座標と y 座標をそれぞれ t の関数として与える T を適当な数 で分割して Δ t とし Δ t 毎の質点の座標を求める これらの座標データをワークシートに書き出す これらのデータを読み込んで グラフを描かせるプログラムは以下のようになる Wth ActveSheet.ChartObjects.Add(50, 0, 00, 00).Chart.ChartType = lxyscatterlesnomarkers.setsourcedata Source:=ActveSheet. _ Rage(Cells(, ), Cells( +, )), _ PlotBy:=lColums.Aes(lCategory).MmumScale = 0.Aes(lCategory).MamumScale =.Aes(lValue).MmumScale = 0.Aes(lValue).MamumScale =.HasLeged = False Ed Wth * 注 : 上記について.Add(50, 0, 00, 00) の括弧内数字は グラフエリアの画面上における左上と右下の座標であり.SetSourceData Source:=ActveSheet.Rage(Cells(,), Cells(+,)) は プロットを行うデータが格納されているセルの範囲である ( 結果を出力したセル範囲に合わせて変更が必要 ) また.Aes(lCategory).MmumScale = 0 等はグラフの, y 軸の表示範囲で 任意に設定された角度 a, 初速度 v 0 に対して 計算結果がすべて表示されるよう変更する必要がある 座標 y 座標 0 0 0.883699.37755.767399.44898.65098 3.486 3.534797 3.67347 0 0 8 6 4 0 5 0 () プログラムリスト フローチャート 質点の軌跡の座標データ 軌跡の図をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 7
演習課題 5 ある学級で行ったテストの結果が下表のように得られている () この表をエクセルで作成し ファイル名を test.ls として保存せよ () test.ls のデータを読み込み 各科目の平均値と標準偏差を計算するプログラムを作成せよ プログラムは test.ls とは異なるファイル名で保存せよ (3) 英語と数学の成績の関係式 ( 次式 ) を最小 乗法により求めよ () 平均値は次式を用いる = = ( = データ数 = 点数 ) 合計は次の反復式より得られる sum = sum + () 標準偏差は次式を用いる σ = = ( ) (3) 最小 乗法 ( 教科書 5.) による解法は以下の通りである 回帰式を y = a + b とおくと f ( a, b) = = { y ( a + b)} が最小となるように a, b を決定する ここでは が英語 y が数学の成績である f = { y ( a + b)} = 0 a = f = { y ( a + b)} = 0 b = 上式をa, b について整理すると σ y a + b = y a = r, b = y a = = = σ = r = a + b = = y = y ( )( y y) σ σ ( 相関係数 ) () プログラムリスト フローチャート 計算結果をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル (3) データファイル (test.ls) テスト結果 出席番号 英語 数学 57 60 93 8 3 66 54 4 4 68 5 5 34 6 77 64 7 69 80 8 7 6 9 88 9 0 9 35 7 6 60 30 3 64 64 4 54 48 5 00 98 6 50 55 7 48 44 8 68 6 9 0 8 0 50 45 57 48 83 88 3 88 90 4 73 70 5 55 56 6 8 7 7 6 65 8 73 76 9 90 84 30 55 8 3 60 5 3 3 55 33 68 74 34 66 56 35 76 68 36 30 50 37 5 5 38 40 44 39 34 38 40 75 7 8
演習問題 6 ある生物種の世代ごとの個体数が, 以下の式で与えられている. { ( ) } (0 4) ( + ) = a ( ) a () ここで, (+) は + 世代目の個体数 ( 子供の数 ), () は 世代目の個体数 ( 親の数 )- () は食料の量を表し,a はその生物の繁殖力を意味している. の第 世代の個体数を 0.5( 億個体 ) として, 繁殖力 a を様々に変えた場合の,30 世代までの個体数の変化をグラフ化して考察せよ. 個体数の変数 を配列変数として扱い,For ~ Net ループを使って計算を行うと良い. すなわち, Dm a, (30) As Double 配列変数 ( 倍精度型 ) の定義 () = 0.5 For = To 30 () =.. 式 () For ~ Net ループ Net なお a は以下の解説に基づき 6 種類以上の値を設定すること 作図については, 作図プログラム例を参考にする ( 例は, データが A3 から G3 の範囲にある場合 ) か あるいはエクセルのグラフウイザードを用いて作成すること 作図プログラム例 Rage("A3:G3").Select Charts.Add ActveChart.ChartType = lxyscatterles ActveChart.SetSourceData Source:=Sheets("Sheet").Rage("A3:G3"), PlotBy _ :=lcolums ActveChart.Locato Where:=lLocatoAsObject, Name:="Sheet" ActveChart.ApplyDataLabels Type:=lDataLabelsShowNoe, LegedKey:=False ActveChart.SeresCollecto().Select ActveChart.Aes(lCategory).Select Wth ActveChart.Aes(lCategory).MmumScaleIsAuto = True.MamumScale = 30.MorUtIsAuto = True.MajorUtIsAuto = True.Crosses = lautomatc.reverseplotorder = False.ScaleType = llear.dsplayut = lnoe Ed Wth 一般に, 式 () はロジスティック写像と呼ばれ, 最も簡単なカオス現象の例として知られている. 繁殖力 a と個体数の関係は以下のようになっている. 入力条件 :0 (),0<a 4 のとき, 0<a< :() は 0 になる. a :() は単調に -/a になる. <a 3 :() は振動しながら -/a になる. 3<a + 6(3.44...):() は周期振動になる. + 6<a 3.57... :() は多重周期振動になる. 3.57...<a 4 :() はカオス現象になる. () プログラムリスト フローチャート 入力データ ( 入力した a の値 ) 30 世代までの個体数の変化の図とその考察をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 9
演習課題 7 常微分方程式の差分解法 ( オイラー法 ) 微生物の増殖速度は以下の式で表現される dx C = X K d X dt μ K c C + ここで X: 微生物濃度 [mg/l] t: 時間 [h] μ: 増殖速度定数 [/h] C: 餌濃度 [mg/l] K c : 定数 [mg/l] K d : 死滅速度定数 [/h] で 右辺第 項が増殖 第 項が死滅を表す また 餌の濃度変化は K Y を定数として C X dc μ K c C + = - dt KY 微生物の初期濃度 (X=0[mg/L]) 増殖速度定数(μ=[/h]) 餌の初期濃度(C=5000[mg/L]) 死滅速度定数 (K d =0.[/h]) 定数 K c =0000[mg/L] K Y =0.9 として 4 時間後の微生物濃度と餌濃度を求めよ 時間刻み (Δt) は 0.5[h] とする C dx μ' = μ とおくと = μ' X K d X となる この式を差分化すると K c + C dt X t+δt X t = μ' X t Δt 従って 時間 t において X t の濃度を有する微生物の t+δt 時間後の濃度 X t+δt は として予測できる X t+ Δt = X t + μ ' X tδt 初期値と各定数をワークシートのセルに記入し 各変数に初期値と定数値を取り込む 微生物と餌の初期濃度および各定数を用いて微生物増殖速度 (μ'x) を求める 3 増殖速度 (μ'x) から餌の消費速度を求める K s = μ' X/ 4 計算した微生物増殖速度を用いて時間ステップ後の微生物濃度を次式で算出する X 0+ Δt = X + μ ' KY ( X K X ) Δt 5 餌の消費速度は 3 から求まるので 時間ステップ後の餌濃度は d C 0 +Δt = C K Δt s から 5 までを ( 出力を行いながら )t が 4 時間になるまで繰り返す 6 エクセルのグラフウイザードを用いて 結果を図にする () プログラムリストと出力結果図をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 30
演習課題 8 ランダムウォークある場所 ( 原点 ) にいる 50 人の酔っ払いが 縦横同じ長さの碁盤状の街路を歩く 酔っ払いは勝手な方向 ( 東西南北 ) に進む * 50 人の酔っ払いが 60 回十字路を進む軌跡を図で示せ * 各方向に進む確率を与える ( 最初は各方向に同じ確率としよう ) 原点 (=0, y=0) の座標を与え (60) 回一様乱数 (0 ~ の間で平均的に一様に出現する変数 ) を生成して進む方向を決め 座標をワークシートに出力する VBA で一様乱数を発生するには 次のステートメントを用いる ( 教科書 p.34) Radomze = Rd ' 乱数系列の初期化 ' に 0~ までの乱数が代入される この手順を c (50) 人分繰り返して座標を ( 行 *c 列 ) 出力した後 c (50) 系列の散布図を作成する 散布図を表示させるプログラムは以下の通りである ActveSheet.ChartObjects.Add(00, 0, 300, 300).Actvate ActveChart.ChartType = lxyscatterlesnomarkers ActveChart.HasLeged = False For j = To c ActveChart.SeresCollecto.NewSeres Wth ActveChart.SeresCollecto(j).XValues = _ Rage(ActveSheet.Cells(, + (j - ) * ), _ ActveSheet.Cells( +, + (j - ) * )).SeresCollecto(j).Values = _ Rage(ActveSheet.Cells(, + (j - ) * ), _ ActveSheet.Cells( +, + (j - ) * )).Aes(lValue).HasMajorGrdles = False Ed Wth Net ここで ( 出力開始行番号 ) c はプログラムに合わせて変更すること () プログラムリスト フローチャート 酔っ払いの軌跡図 考察 ( 何度か計算を繰り返し 結果の違いを観察する 上下左右に進む確率を変更して実験を行ってみるなど ) をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 3
演習課題 9 ガウス ザイデル法による連立 次方程式の解法ガウス ザイデル法により 以下の 3 元連立 次方程式を解くプログラムを作りなさい 繰り返しは 0 回とし 初期値によって解がどのように収束していくのか確認すること そして 初期値による収束の様子の違いをレポートにまとめよ 0 + 3 = 6 0 33 = 7 () 3 + 5 + 0 = 99 3 [ 解説 ] ガウス ザイデル (Gauss-Sedel) 法は 連立一次方程式の解法のひとつで 適当な初期値から解を収束させていく方法である 教科書にある 掃き出し法 ( 教科書 第 4 章 ) のような代数的な解法ではなく またいつでも収束するわけではないが 演算が速く 多変数の場合にも向いている 以下の連立方程式を例にアルゴリズムを説明する 3 + + 3 = 8 + 6 + 3 3 = 3 + + 4 3 = 7 まず これらを = ~, = ~, 3 = ~ という形に変形する = ( 8 - - 3 ) / 3 = (3 - - 3 3 ) / 6 3 = (7 - - ) / 4 ここで適当な解 ( 値は何でも良い ) を設定する 仮に (,, 3 ) = (,,) とすると = (8--)/3 = となり さらに第二式より = (3-4-3)/6 =.666 3 = (7--5.333)/4 =.4667 となる これを繰り返すと 3 回目 0.97.30093.85648 3 回目 0.94753.0895.96849 4 回目 0.980753.07.99373 5 回目 0.99470.00490.9988 6 回目 0.99874.00098.9998 7 回目 0.99973.0008.99998 8 回目 0.999948.00003 3 9 回目 0.99999 3 0 回目 0.999999 3 回目 3 となり 解が求まる この課題では 0 回繰り返すことになっているが 通常のプログラムでは 許容誤差を前もって与え 繰り返しによる解の変化がそれよりも小さければ収束したとして終了する 3
以下の連立一次方程式を考える a +a + +a = b a +a + +a = b a +a + +a = b ガウス ザイデル法の漸化式は次のように与えられる ( k+ ) ( k+ ) ( k ) = b aj j aj j / a j= j= + (k+) ここで a j および a b : 繰り返し回数 k+ 回目における ( =,,3) の解 : 式 () 左辺の係数項 : 右辺定数項 である ( 参考 ) ガウス ザイデル法が収束するための十分条件は 以下の式で与えられる a a j j= ガウス ザイデル法の計算ステップ : パラメータとデータの入力 : 方程式の数 ( = 3) 反復回数 kc (kc = 0) 方程式の係数項を入力する 方程式の係数項は配列で取り扱う : 番目の方程式を対角要素 a (a 0) で割る a j =a j / a ( =~;j = ~+) 3: 初期近似解 を任意の数 (0) に設定する 4: 上記の漸化式を使って 逐次近似解 k+ を kc 回反復計算する 5: 解を出力する = k+ ( =~; = 3) 考察は, 初期値によってどの程度の繰り返し回数で理論解に達するか検討すること 初期値は 3 種類以上について比較を行い 図を用いて説明を加えてもよい () プログラムリスト フローチャート 出力結果 初期値と収束の様子との関係についての考察をまとめた Word ファイル () 作成したプログラムを含む Ecel ファイル 33