[size=150][b]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][br][/b]前回まで、微分方程式の基本となる1階線形の解法を扱ってきた。[br]今回は、同じ線形だけど、2階の場合について調べてみよう。[br][b][br][size=150][b]＜定数係数の場合＞[/b][/size][br][br]2階線形微分方程式の基本は定数係数だ。[br][br][/b][/size]２階のa、bが定数の線形微分方程式[b]y"+ay'+by=0[/b]は2次方程式の解法に似ている。[br]等比数列の漸化式の単純化や行列の固有値を求めるのにも登場した[br]あの固有方程式、特性方程式が登場する。[br][br][b][color=#0000ff]固有方程式（特性方程式）[/color][/b]t[sup]2[/sup]+at+b=0の解がt=p,qのとき、[br]一般解は固有方程式の解の種別によって３タイプある。[br]・p,qが異なる実数なら[b]y=Ae[sup]px[/sup]+Be[/b][sup][b]qx[/b][br][/sup]・p=qの重複実数なら(p=-a/2)で[b]y=(A x +B)e[/b][sup][b]px[/b][br][/sup]・複素数解p=m+ni, q=m-niなら、y=Ae[sup]px[/sup]+Be[sup]qx[/sup]＝[b]e[sup]mx[/sup](C cosnx + D sinnx)[br][/b]　m=-a/2, n =√D/2[br][color=#0000ff][b]（理由)[/b][/color][br]y'+ky=0の解は、dy/dx=-kyを-1/ky dy=1 dxで変数分離すると、-1/k∫1/y dy=∫1 dxとなる。[br]-1/k log|y|=x＋cから、 log|y|=-kx +d Cを定数とする一般解はy=Ce[sup]-kx[/sup]と指数関数になる。[br]だから、[b]y"+ay'+by=0 の１つの解をy=e[sup]tx[/sup]とおいてみる。[br]y'=t e[sup]tx[/sup] ,y'' = t[sup]2[/sup] e[sup]tx[/sup][/b]となる。[br]これをもとの方程式に代入すると、 t[sup]2[/sup] e[sup]tx[/sup]+a *t e[sup]tx[/sup] +b*e[sup]tx[/sup] =(t[sup]2[/sup] + at +b)e[sup]tx[/sup] となる。[br]だから、tが特性方程式(t[sup]2[/sup] + at +b)=0　の解なら、y=[b]e[sup]tx[/sup][/b]は微分方程式の解だといえるね。[br][br]・特性方程式が異なる実数解p,qを持つなら、２関数y1=e[sup]px[/sup], y2=e[sup]qx[/sup]が微分方程式の解だから、[br]　微積分の演算はベクトルのように線形だから、[br]　この2関数を基底とした線形結合y=Ay1+By2が一般解となるね。[br]　つまり、[color=#0000ff][b]y=Ae[sup]px[/sup]+Be[/b][sup][b]qx[/b][/sup][/color][br]・特性方程式が重複実数解をもつなら判別式＝０で、t=-a/2が重複解。[br]　y1=e[sup]-ax/2[/sup]、もう一つの解はy2=ｘ*y1とするとy1に比例しない独立な基底にできる。[br] だから、[b][color=#0000ff]y=Ay1+By2=Ae[sup]-ax/2[/sup] +Bxe[sup]-ax/2[/sup]=(A+Bx)e[/color][/b][sup][b][color=#0000ff]-ax/2[/color][/b][br][/sup]・特性方程式が複素数解p=m+ni, q=m-niなら、[br] P=e[sup]px[/sup]= e[sup]mx+inx[/sup]=e[sup]mx [/sup]e[sup]inx[/sup]=e[sup]mx [/sup](cos nx + i sin nx) [br] Q=e[sup]qx[/sup]= e[sup]mx-inx[/sup]=e[sup]mx [/sup]e[sup]i(-n)x[/sup]=e[sup]mx [/sup](cosnx - i sin nx)[br] (P+Q)/2=e[sup]mx [/sup](cos nx)=y1 , [br] (Q-P)/2i=e[sup]mx [/sup](sin nx)=y2　とおくと、y2/y1=tan nx で、[br] 一定ではないので独立基底とみることができる。[br]　[b]一般解は[color=#0000ff]y=AP+BQ=Cy1+Dy2=e[sup]mx [/sup](C cos nx+D sin nx)[/color][/b][br][color=#0000ff]（例）[br][/color]「y" -3y’+2y=0の一般解」は？[br]固有方程式(t-1)(t-2)=0の解t=[b]1,2[/b]から、y=Ae[sup][b]1x[/b][/sup]+Be[sup][b]2x[/b][/sup][br][color=#0000ff]（例）[br][/color]「y" +2y’+y=0の一般解」は？[br]固有方程式(t+1)[sup]2[/sup]=0の解t=[b]-1[/b]から、y=(Ax+B)e[sup][b]-1x[/b][/sup][br][color=#0000ff]（例）[/color][br]「y" -4y’+13y=0の一般解」は？[br]固有方程式の解t=[b]2[/b]±[b]3[/b]iから、y=e[sup][b]2x[/b][/sup](A cos[b]3x[/b] + B sin[b]3x[/b])[br]

[color=#9900ff][u][b][size=150]質問：定数係数の2階線形微分方程式を解析的に解くコードはどうかけばよいでしょう。[br][/size][/b][/u][/color][br][b][size=150]＜Sympy＞[br][/size][/b]sympyで、y(n)(x)はDerivative(y(x), x, n)とかけます。[br]だから、[br]y''を変数ddyにDerivative(y(x),x,2)として入れて、[br]y'を変数dyにDerivative(y(x), x)として入れておくと、[br]微分方程式Eq(左辺, 0 )の左辺の中に使えます。[br]「y" -3y’+2y=0」y=Ae[sup][b]1x[/b][/sup]+Be[sup][b]2x[/b][/sup][br]「y" +2y’+y=0」y=(Ax+B)e[sup][b]-1x[/b][/sup][br]「y" -4y’+13y=0」y=e[sup][b]2x[/b][/sup](A cos[b]3x[/b] + B sin[b]3x[/b])[br]を確認しよう。[br][br]#[IN]Python[br]#===========================[br][color=#0000ff][b]from sympy.abc import *[br]from sympy import *[br]y = Function('y')[br]ddy=Derivative(y(x), x, 2)[br]dy=Derivative(y(x), x)[/b][/color]

[b][size=150]＜geogebra＞[br][/size][/b]geogebraでも常微分方程式が解けます。[br]それは、[url=https://geogebra.github.io/docs/manual/en/commands/SolveODE/]SolveODE[/url]です。[br][br][b][color=#0000ff]SolveODE(b(x), c(x), f(x), Start x, Start y, Start y', End x, Step)[br]Solves second order ODE y′′+b(x)y′+c(x)y=f(x).[br][/color][/b]とありますね。[br]「y" -3y’+2y=0」y=Ae[sup][b]1x[/b][/sup]+Be[sup][b]2x[/b][/sup][br]「y" +2y’+y=0」y=(Ax+B)e[sup][b]-1x[/b][/sup][br]「y" -4y’+13y=0」y=e[sup][b]2x[/b][/sup](A cos[b]3x[/b] + B sin[b]3x[/b])[br]を確認しよう。[br][br][b]SolveODEは、2階の場合は数値計算をして線をつないでいる[/b]ようです。[br]1階ならば、設定がかなり少なくなるのは、解析的に解いて定数を代入するしくみなのでしょうね。[br]作った例が冒頭の３つのアプレットです。[br][b]SolveODE(b(x), c(x), f(x), Start x, Start y, Start y', End x, Step)[br][/b]SolveODE(-3,2,0,-3,0.1,0.01,1,0.1)[br]SolveODE(2,1,0,-1,-4,8,4,0.1)[br]SolveODE(-4,13,0,-5,0,0.001,3,0.1)[br]と設定しておきました。[br]一般解の式のA,Bをうまくさがすと、ピッタリグラフが重なりますね。[br][br]同じgeogebraでもCASのSolveODEでは、解析に式を返します。[br]これでやってみましょう。[br][b][color=#0000ff]Sympyと同様に解析的に数式を返してくれます。[br]CASはすばらしい。[/color][/b]

[b][size=150]＜ばねの減衰しない振動＞[br][/size][/b][br]ばねの長さの変化をしらべよう。[br][b]フックの法則F=－kx[/b](ばねが自然長からｘだけ引っ張るときに生まれる力Fはｘに比例する）[br][b]ニュートンの運動方程式F＝ｍa[/b](等速度運動か静止を変更する力は物の質量と加速度の積に等しい）[br]この２つから、微分方程式を作りたい。つまり、空気抵抗などの邪魔がないときを考えてみよう。[br][br]変位ｘは時間tの関数で、加速度a=x''だから、mx''=-kx[br]つまり、x''+k/m x=0という2階線形微分方程式ができたね。[br]k/m=ω[sup]2[/sup]とおくと。[b][color=#0000ff]x''+ω[sup]2[/sup] x=0[/color][/b]です。[br]初期値は、ばねの自然長さがx0、t=0の速度はv=0としましょう。[br][br]・特性方程式t[sup]2[/sup]+ ω[sup]2[/sup]=0の解は複素数解p,q=0±i ω[br]一般解はx=e[sup][b]0x[/b][/sup](A cos[b]ωt[/b] + B sin [b]ωt[/b])[br][b][color=#0000ff]x= A cosωt + B sin ωt[br][/color][/b][br]・確認してみると、[br]v=x'=[b]ω[/b] (A -sin[b]ωt[/b] + B cos [b]ωt)[br][/b]a=x''=-[b]ω[/b] [sup]2[/sup](A cos[b]ωt[/b] + B sin [b]ωt)=-[b]ω[/b] [sup]2[/sup]x[br][/b][color=#0000ff][b]x''=-ω [sup]2[/sup]x[/b][/color]が出ました。[br]これで、一般解がもとの微分方程式を満足させていることがわかりました。[br][br]・初期値はまだ使ってませんでした。[br]一般解x= A cos[b]ωt[/b] + B sin [b]ωt[/b]で、[br]t=0と、x=x0から、x0=A cos[b]0[/b] + B sin [b]0=A[/b]速度v=[b]ω[/b] (A -sin[b]ωt[/b] + B cos [b]ωt)[/b]で、[br]t=0と、v=0と、A=x0から、0=[b]ω[/b] (x0 -sin[b]0[/b] + B cos 0)=ω B 、ω≠0なので、B=0[br]ということは、x= A cos[b]ωt[/b] + B sin [b]ωt[/b]= x0 cosωt[br][b][color=#0000ff]特殊解はx=x0 cosωt[br]このとき、v=x'=-ωx0 sinωt[br][/color][/b]ωt=2πを解くと、t=2π/ωですから、[br]この特殊解は、振幅x0で周期T=2π/ωの振動を表していますね。[br]振動数（周波数）は周期Tの逆数で、1秒間に1/T=ω/2π振動します。[br]単振動とか調和振動というものです。[br][br]・さらに、特殊解を分析してみよう。[br][b]v[sup]2[/sup][/b]=(-ωx0 sinωt)[sup]2[/sup]=(ωx0)[sup]2[/sup](sinωt)[sup]2[/sup]=(ωx0)[sup]2[/sup](1- (cosωt)[sup]2[/sup])= ω[sup]2[/sup](x0[sup]2[/sup]-(x0cosωt)[sup]2[/sup])=[b]ω[sup]2[/sup](x0[sup]2[/sup]-x[sup]2[/sup])[/b]だから、[br]v[sup]2[/sup]=(k/m)(x0[sup]2[/sup]-x[sup]2[/sup])[br]両辺にm/2をかけて、[b]1/2mv[sup]2[/sup]= 1/2kx0[sup]2[/sup]-1/2kx[/b][sup][b]2[/b][br][/sup]運動エネルギーがばねの弾性エネルギーの減少分に等しいことを言っています。[br]つまり、ばねについての力学的エネルギーが運動エネルギーと弾性エネルギーの和であり、保存されるということを言っているのですね。

[b][size=150]＜ばねの減衰振動＞[br][/size][/b][br]空気抵抗を考慮するとどうなるだろうか。[br]空気抵抗は、ばねの先端のおもりの速さに比例して、変位と逆向きに働くから-cx'とおけるね。[br]mx''=-kx-cx'(m、k、ｃも正）Cは抵抗の係数になっている。[br]つまり、[b]mx''+cx'+kx=0[/b]という2階線形微分方程式ができる。[br]初期値は、ばねの自然長さがx0、t=0の速度はv=0としましょう。[br][br]・特性方程式mt[sup]2[/sup]+ct +k=0の解は場合分けが必要になる。解がt=p,qのとき、[br]一般解は固有方程式の解の種別によって３タイプある。[br]判別式はD=c[sup]2[/sup]-4mk[br]①型・c[sup]2[/sup]-4mkが正で、p,qが異なる実数ならx[b]=Ae[sup]pt[/sup]+Be[/b][sup][b]qt[/b][br][/sup]②型・c[sup]2[/sup]-4mk=０で、p=qの重複実数なら(p=-c/2m)で[b]y=(A t +B)e[/b][sup][b]pt[/b][br][/sup]③型・c[sup]2[/sup]-4mkが負で、複素数解p=a+bi, q=a-biなら、y=Ae[sup]pt[/sup]+Be[sup]qt[/sup]＝[b]e[sup]at[/sup](C cosbt + D sinbt)[br][/b]　a=-c/2m, b =√D/2[br][color=#0000ff][br]動きの特徴[br][/color]①型・c[sup]2[/sup]-4mkが正で、Cがどでかい。[br]m、c、kが正の2次関数f(t)=mt[sup]2[/sup]+ct+kとf=0との交点がt=p,qになるとしよう。[br]f(0)=kは正。軸t=(p+q)/2= -c/2mは負。だから、p,qともに負となるね。[br]一般解は、x=Ae[sup]pt[/sup]+Be[sup]qt[br][/sup]の指数部分pt, qtはｔが正の時間でつねに負となる。ということは、ｔの増大により[br]exp(pt), exp(qt)は無限に０に近づく。[br]振動することなく、空気抵抗で押しつぶされる。[br][br]②型・c[sup]2[/sup]-4mk=０で、Cがそこそこ大きい。[br]m、c、kが正の2次関数f(t)=mt2+ct+kとf=0との交点がt=p=qになるとしよう。[br]f(0)=kは正。軸t=(p+q)/2= -c/2mは負。だから、t=p=qは負となるね。[br]一般解はy=(A t +B)e[sup]pt [/sup]①型と同様に振動することなく縮む。[br][br]③型・c[sup]2[/sup]-4mkが負で、ｃは小さいので単振動に近いが、ｃは０ではないから減衰がある。[br]複素数解p=a+bi, q=a-biなら、y=Ae[sup]pt[/sup]+Be[sup]qt[/sup]＝e[sup]at[/sup](C cosbt + D sinbt)[br]a=-c/2m, b =√-D/2[br]一般解y=exp(at)( C cosb[b]t[/b] + D sin b[b]t)[/b]で、t=0と、y=x0から、[br]x0=exp(0)(Ccos0 +Dcos0)=C。[br]速度v=y'=(e[sup]at[/sup])'(C cosbt + D sinbt)+(e[sup]at[/sup])(C cosbt + D sinbt)'[br]=a(e[sup]at[/sup])(C cosbt + D sinbt)+b(e[sup]at[/sup])(-C sinbt + D cosbt)[br]=exp(at){ (aC+bD)cosbt +(aD-bC)sinbt} [br] t=0と、v=0と、C=x0から、[br]0=exp(0){ (ax0+bD)cos0 +(aD-bx0)sin0} =(ax0+bD)なので、D=-x0 * a/b。[br]ということは、特殊解はx= e[sup]at[/sup](x0 cosbt -x0 *a/b sinbt)[br][color=#0000ff][b](例)[/b][/color][br]初期値y0=0.15(m), 初速v0=0[br]質量m=89/9.8=9.082(kg),[br]復元定数k =890(kg/s[sup]2[/sup]) ,[br]抵抗係数c=0から200としたときの[br]ばねの動きを予想しよう。[br][br][b]・c=0のときは単振動（調和振動）[br][/b]k/m=890/(89/9.8)=98=ω[sup]2[/sup][br]ω=9.899[br]ω/2π＝9.899/6.28=1.57555[Hz] 1秒の振動数[br][b][color=#0000ff]y=y0 cosωt=0.15 cos(9.899 x) が特殊解。[br][/color][/b][br]・c=200以下で減衰振動する範囲は？[br]f(t)=mt[sup]2[/sup]+ct +k=9.082 t[sup]2[/sup]+ct +890=0が複素数解をもつとき。[br]D=c[sup]2[/sup]- 4*m*kが負, D' =-D[br]s=pow(4*k*m,1/2)[br]s=179.81078944268054[br][b]cが179.81078944268054以下のとき減衰振動になる。[br][/b][br]c=100のとき、[br]a=-c/2m=-5.50539[br]b=√-D/2=pow(-c**2+4*k*m,1/2)/2m=8.22719[br]C=x0=0.15[br]D=-x0*a/b=0.10038[br]特殊解は[br][b]y= e[sup]at[/sup](x0 cosbt -a x0 sinbt)[br][/b][b]=exp(-5.50539 x) (0.15 cos( 8.22719 x) + 0.10038 sin( 8.22719 x))[br][/b]