振動する梁のたわみ関数その3 - Maximaでこうぞうりきがく

以前，片持ち梁の固有振動数を求める際に一次モードの梁のたわみ関数を与えましたが
簡単な調和振動関数を前提としたもので精度はそれなりでした。
今回はもう少しちゃんと計算をしたら？というお話です

%i1にてパッケージ"newton1"をロードします（画面出力は省略）
newton1はニュートン法を利用するためのパッケージです（ニュートン法その1を参照ください）
梁のたわみ関数y(x)を任意の調和振動関数(%o2)で仮定します
振動する梁のたわみ関数その1に比べて係数CとDの項が増えてます

【片持ち梁】

片持ち梁の境界条件を指定します（係数が2つ増えたので境界条件も2つ増えてます）

x = 0 の位置の変位が 0 (%o3)
x = 0 の位置のたわみ角が 0 (%o4)
x = l の位置の曲率が 0 (%o5)
x = l の位置の3階微分が 0 (%o6)

上4式を未定係数A～Dに対しての連立方程式として，係数マトリックスを%o7式に示します
係数マトリックスの特異性より%o8式を得ます
この方程式を ω*l について解きたいのですが残念ながら陽には解けません
仕方ないのでニュートン法を使って求めた（1次の振動モードに対応する）数値解を%o9式に示します
1次以外の数値解を求めたい場合は初期値(x0 = 2)を別の値へ修正してください
余談ですがこの ω*l を固有値と呼びます

係数マトリックスと未定係数ベクトルの積を%o9式でまとめた結果が%o11式となります
上式の2～4成分を連立方程式としてB～Dについて解いた結果を%o12式に示します
この解で%o2式を書き直したものを%o13式に示します（Aは任意の振幅）