X''=(Y-CX'-BX)/A パラメータの導出

質問

X',X''はそれぞれXの一回微分，二回微分を示します。

X＝0.01，Y＝0.39
X＝0.008，Y＝0.37
X＝0.09，Y＝0.36

という値がある時，パラメータ，A,B,Cを求めたいのですが，どのように求めるのでしょうか？
よろしくお願いします。

回答No.11

info22

>x=C1exp(α1y)+C2exp(α2y)　　・・・(1)
α1,α2が虚数になる時も同次微分方程式の解は
形式上はこの式になりますが、
xは実数なので(1)の虚数を含んだ表現は適当でないので、実数だけで表した解の式を使うべきです。
虚数を最初から含んだ式を使えば、最終的に虚数を含まない式にオイラーの公式を使って戻す必要があり、未習熟な方は混乱するだけです。
(1)の代わりに次の(1)'の式を使えば実数だけのyの式になり、以降実数だけの計算になります。

x=exp(-cx/(2a))[C1'*cos(√{4*a*b-c^2))*y/(2*a)}
+C1'*sin(√{4*a*b-c^2))*y/(2*a)}] …(1)'

>#5のx(0)=x0(初期位置),dx(0)/dy=v0(初期速度)の初期条件より

確認ですが、元の微分方程式のYに相当するy、Xに相当するxと初期値x0,
v0との関係はどうなっていますか？
以下の３つのケースのいずれですか？

(i) x(y),yは独立変数（yは時間）
　これまでの回答はこのケースで話が進んでいます。
(ii) x(t),y(t),tは時間
(iii) x(y),y(t),tは時間
　　(ii)と(iii)のケースだとy(t)を与えないと
　　非同次方程式の特解が求まらないので
　　もとの微分方程式自体が解けないことになります。

(i)のケースを前提に回答をします。

この初期条件は
>x=C1exp(α1y)+C2exp(α2y)　　・・・(1)

この式でなく
x=[exp{-cy/(2a)}]〔C1'cos[{√(4ab-c^2)}y/(2a)]+C2'sin[{√(4ab-c^2)}y/(2a)]〕
+(yb-c)/(b^2) …(1)"
に適用しないと意味がないですね。

>x0=C1+C2 ・・・(2)
>v0=-α1C1-α2C2　　・・・(3)
>(2),(3)より
>C1=v0/-α1+α2
>C1=a*v0/-i√(c^2-4*a*b)
>C2=a*v0/i√(c^2-4*a*b)
>と求めてしまいました。

以上は次のように書き換わります。

y=y0=0
x0=x(y0)=C1'-(c/b^2)
v0=dx/dy(y0)=-{c/(2a)}C1'+[{√(4ab-c^2)}/(2a)]C2'+(1/b)

C1'=x0+(c/b^2)
C2'=[v0-(1/b)+{c(x0+(c/b^2))/(2a)}]/[{√(4ab-c^2)}/(2a)}
={2av0+(2a/b)+cx0+(c/b)^2)}/√(4ab-c^2)

これを上の(1)"に代入すればいいですね。

以下は本来、虚数でない式に形式上虚数を使って表しただけ式ですから、
(1)"と同じ式です。わざわざ逆行して虚数表現式のC1,C2に戻すことは全くナンセンスです。

>X=C2*exp[-{C-i√(4AB-C^2)}Y/(2A)]+C1*exp[-{C+i√(4AB-C^2)}Y/(2A)]+(YB-C)/(B^2)
>のα1，α2とも対応していると思ったのですが…
>任意定数を(C1、C2)から(C1',C2')に置き換えたとありますが、
>それぞれC1=　C2=　と表すとどのようになるのでしょうか？
>虚数に関しての理解が浅いため、C1',C2'と結びつけた、境界条件の与え
>方で虚数計算を出なくしたいです。

上で求めたC1',C2'を(1)"に代入した式が
c^2<4abのケースにおける
元の非同次微分方程式の（虚数を含まない）解xの式になります。

補足

ありがとうございます。

(i)のケースです。
(1)"に、C1',C2'を代入してフィッティングを行ってみたところ、データによって-になってしまった場合と
すべて+で出た場合とがあったので、今後#9を参考にデータをより正確に計測してフィッティングを行っていきたいと
思います。

gnuplotのパラメータの初期値の設定ですが、式から判断できない場合、
すべて1で行ってもよいのでしょうか？
その場合、一度目のフィッティングで求められた値を初期値に用いて再度フィッティングを行うというような必要はあるのでしょうか？

ご回答よろしくお願いします。

回答No.10

info22

#1～#9です。
>への求め方は、オイラーの公式を用いているのでしょうか？
そうです。
そして、任意定数を(C1、C2)から(C1',C2')に置き換えただけです。
C^2<4ABの場合の微分方程式の解を最初から
X= =[exp{-CY/(2A)}]〔C1'cos[{√(4AB-C^2)}Y/(2A)]+C2'sin[{√(4AB-C^2)}Y/(2A)]〕
+(YB-C)/(B^2)
の形で書いてある教科書や参考書が殆どなので、まわりくどくオイラーの公式を使う必要が無いかもしれません。A#2では機械的にXを求めると虚数が指数部に出てくる→その場合はオイラーの公式でsin,cosの式に直せば虚数単位iが解の式からなくなることを示しただけですので、余り気にしないで下さい。やりたければ、自分でゆっくり時間をかけてやってみてください。

>また、Q#6の式よりC^2<4ABの場合の境界条件を求める過程なのですが、
>α1=(-c+i√(4*a*b-c^2))*x/(2*a)
>α2=(-c-i√(4*a*b-c^2))*x/(2*a)
なぜ本来、実数であるはずの境界条件として、
α1,α2を虚数として定義するのですか？

そのため
>と求めました。しかし、i(虚数)が含まれてしまうのですが、上と同じようにiを取り除く計算方法はありますでしょうか？

ということになります。

最初から
X=[exp{-CY/(2A)}]〔C1'cos[{√(4AB-C^2)}Y/(2A)]+C2'sin[{√(4AB-C^2)}Y/(2A)]〕+(YB-C)/(B^2)
の解の形に即した、（C1,C2は使わない。C1,C2を使うから虚数単位iが式中に出てきて混乱します。もちろん理解が深ければ虚数単位iが出てきてもC1,C2をC1',C2'置き換える操作をすれば実数に戻りますが...）、C1',C2'と結びつけた、境界条件の与え方をすれば何もあなたの言う境界条件から虚数を取り除いて、C1',C2'を導く問題は発生しません。
そうすれば、虚数計算が最初から出てこなくなります。

補足

ご回答ありがとうございます。

>実数であるはずの境界条件として、α1,α2を虚数として定義するのですか？

#5のx(0)=x0(初期位置),dx(0)/dt=v0(初期速度)の初期条件より

y=C1exp(α1t)+C2exp(α2t)　　・・・(1)
x0=C1+C2 ・・・(2)
v0=-α1C1-α2C2　　・・・(3)

(2),(3)より
C1=v0/-α1+α2
C2=v0/α1-α2

として

C1=a*v0/-i√(c^2-4*a*b)
C2=a*v0/i√(c^2-4*a*b)

と求めてしまいました。また、

X=C2*exp[-{C-i√(4AB-C^2)}Y/(2A)]+C1*exp[-{C+i√(4AB-C^2)}Y/(2A)]+(YB-C)/(B^2)

のα1，α2とも対応していると思ったのですが…

任意定数を(C1、C2)から(C1',C2')に置き換えたとありますが、
それぞれC1=　C2=　と表すとどのようになるのでしょうか？

虚数に関しての理解が浅いため、C1',C2'と結びつけた、境界条件の与え方で虚数計算を出なくしたいです。

質問ばかりで申し訳ありません。
よろしくお願いします。

回答No.9

info22

#1～#8です。

>A#2の
>C^2<4ABの場合
>でフィッティングを行おうと思うのですが、C1'cos、C2'sinのそれぞれ
>C1'、C2'とは、境界条件をそのまま適用してもよいのでしょうか？

フィッティング前のC^2<4ABの場合の微分方程式の解の式に
そのまま境界条件を代入して決めれば良いですね。

>また、桁落ちに関してなのですが、これ以上の精度でデータを取れない
>場合、取得したデータ上の問題ということで対処のしようがないので
>しょうか？

◆エクセルで計算する場合エクセルの中の計算精度は14～15桁あると思いますが、表示する小数以下の桁数を最大限にしておかないと、
0.000000000000123
といった表示になってそのままa,b,c等の値として使うと計算の有効桁数が3桁になってしまったりします。また計算式の各項に対して有効桁数の少ない項の絶対値が大きい場合は見かけ上の計算値の桁数が14桁でも、有効桁数は有効桁数が絶対値の大きい項に左右されてしまいます。
なので、係数定数a,b,cの全てについて、先頭の0を除いた数値部分の桁数
を同じ桁数分（例えば10桁～14桁）正しく導出しておかないといけません。また絶対値が殆ど同じで符号反対の項を加える計算を行う場合も有効桁数が落ちますので、下のa,b,cの有効桁数は大きい方が望ましいですね。
エクセルが表計算ソフトで一通りの数式計算機能はありますが、専用の数式処理ソフトではありませんので、計算精度は専用ソフトより劣ります。
なのでより有効桁数の計算を行うには数式処理ソフト（gnuplot,R、その他）を使った方が良いですね。

エクセルは万能のソフトですが個々の計算機能は専用ソフトの機能には劣り、計算精度も余り無いので、計算結果に対してはチャンと検証しておきたいですね。

◆エクセルで計算の桁落ちがする場合は

計算精度の取れる数式処理ソフトを使用すると良いですね。
数式処理ソフトは有効桁数を設定して何桁の計算精度の計算も可能ですので数式処理ソフトで有効桁数を20桁以上に設定してフィッティングをしより精度の高いa,b,c,dの値を計算して、必要な有効桁数だけ使えばいいかと思います。

A#7の補足に書かれたa,b,cは
>C^2<4ABの場合
のフィッティング関数を使って求めたものではありませんので
改めて
C^2<4ABの場合
の微分方程式の解のフィッティング関数を用いてフィッティングを行ってください。そして求まったa,b,cに対して場合の条件
c^2<4ab
を満たしていれば正しくフィッティングが行えた事になる行えた事になります。
しかし、この場合も、(xi,yi)のデータの全部が「C^2<4ABの場合」を満たすデータで他の場合のデータが殆ど混入していないことが条件です。

これを確認するにはデータの組を２～３つ位のグループに分割して、個々のグループのデータに付いて同じようにフィッティングを行って、同じ場合分けの同じフィッティング関数が該当し、a,b,cもほぼ同じになれば、分割しないもとのデータのフィッティング関数が適切に求められていることが裏付けられます。

補足

ご回答ありがとうございます。

C^2<4ABの場合の
X=C2*exp[-{C-i√(4AB-C^2)}Y/(2A)]+C1*exp[-{C+i√(4AB-C^2)}Y/(2A)]+(YB-C)/(B^2)

から、

=[exp{-CY/(2A)}]〔C1'cos[{√(4AB-C^2)}Y/(2A)]+C2'sin[{√(4AB-C^2)}Y/(2A)]〕+(YB-C)/(B^2)

への求め方は、オイラーの公式を用いているのでしょうか？
詳しいご説明をお願いしたいです。

また、Q#6の式よりC^2<4ABの場合の境界条件を求める過程なのですが、

α1=(-c+i√(4*a*b-c^2))*x/(2*a)
α2=(-c-i√(4*a*b-c^2))*x/(2*a)

これを次の式に代入して

C1=v0/-α1+α2
C2=v0/α1-α2

C1=a*v0/-i√(c^2-4*a*b)
C2=a*v0/i√(c^2-4*a*b)
（C1,C2は境界条件）

と求めました。しかし、i(虚数)が含まれてしまうのですが、上と同じようにiを取り除く計算方法はありますでしょうか？

桁落ちにつきましても、丁寧にご回答していただきましてありがとうございます。
確かに、他の場合のデータが混入していないという条件は大切ですね。桁落ちと他の場合のデータの混入に気をつけてフィッティングを行いたいと思います。

ご回答よろしくお願いします。

回答No.8

info22

#1～#7です。

A#7の補足の質問について
>フィッティングの結果、
>a = 0.000508
>b = 0.265508
>c = -0.0185804
>となり、ｃが負の数となってしまうのが気になります。

どのフィッテング関数を使い、その条件に合う(Xi,Yi)を使うかが問題です。
なので、A#2に書いたa(=A),b(=B),c(=C)の場合分けに対応した微分方程式の解をフィッティング関数に選ばないといけません。
その関数形はフッテッング用データをプロットしたとき、そのデータセット{(Xi,Yi)｜i=1,2,...N}が
微分方程式のどの解のグラフの形に近いかによって、フィッテング関数の形を選んでやらないといけません。
それを間違えると意味の無い結果が出てくるだけです。
従って、データセット{(Xi,Yi)｜i=1,2,...N}をプロットして、そのデータの散布状況がA#2のどの微分方程式の解の関数のタイプになるかをまず判断し、その解関数をフィッテッング関数としてフィッティングして下さい。
つまり、あなたが使ったフィッテング関数とそのときの(a,b,c)の条件が一致していないと、使ったフィッテング関数がデータセットを代表していないことになります。

あなたがやったフィッテングの結果のa,b,cから
A#2の場合分けの条件を調べて見ると
c^2-4ab= -0.000194...＜0
となります。
これは A#2の場合分けの
>C^2<4ABの場合(a=A,b=B,c=C,C1=d=D,C2=e=Eに対応)
になります。
使ったフィッティング関数がこの場合のものを使っていれば正しいフィッティング結果が得られたことになります。
そうでなければ、正しい結果とは言えません。
求めたa,b,cを使って微分方程式を解くと
x(y) = exp(18.28740158*y)*sin(13.71949296*y)*C2
　　　+exp(18.28740158*y)*cos(13.71949296*y)*C1
　　　+.2635666685+3.766364855*y
となります。
C1とC2は微分方程式の境界条件から求まる任意定数です。

もう一つ、気になることは、
フィッテング計算での有効桁数で
>a = 0.000508
>b = 0.265508
>c = -0.0185804
aの有効数字が3桁しかないことで、
場合分けの判別
D＝c^2-4ab= -0.000194...＜0
の値（絶対値の値）が小さく（計算精度も3桁しかないaを使っているため）
Dの計算が桁落ちして判別が正確にできているかが気になります。

補足

ありがとうございます。
このような意味での場合分けだったのですね！
ようやく理解できました。

A#2の
C^2<4ABの場合
でフィッティングを行おうと思うのですが、C1'cos、C2'sinのそれぞれC1'、C2'とは、境界条件をそのまま適用してもよいのでしょうか？

また、桁落ちに関してなのですが、これ以上の精度でデータを取れない場合、取得したデータ上の問題ということで対処のしようがないのでしょうか？

ご回答よろしくお願いします。

回答No.7

info22

#1～#6です。

>d(D),e(D)の式のxは#1のxの値に対応しているのですよね？
？
d(D),e(E)と書いたはずですが？

dとDとC1, eとDとC2
は同じ意味（内容）の任意定数で、使った式で文字を混同しないように、あるいはgnuplotで小文字を使うこともあり、質問者さんの式でも別の文字で使われていたので、同じ意味（内容）の定数を単に別の文字を使っているだけということを指摘しただけに過ぎません。

d,eは微分方程式の初期条件（境界条件）から決まる任意定数ですので、フィッティングの対象となる文字係数は,A,B,C(あるいはa,b,c)の３つですね。
なのでA#4で書いた文字係数a,b,c,d,eの内のd,eは境界条件から決める任意定数なので、d,eを除いた文字係数a,b,cの３つでフィッティングをして下さい。

>エクセル(下記参照)でもできるようなのでやってみたいと思います。
当方はエクセルではやったことがないですが、やる場合は、計算精度の桁落ちが起こす過去の質問が見られますので、計算の有効桁数を十分確保して行って下さい。

補足

ありがとうございます。

エクセルでもgnuplotでも求めることができました！

フィッティングの結果、
a = 0.000508
b = 0.265508
c = -0.0185804

となり、ｃが負の数となってしまうのが気になります。（a,b,cはスカラー量のため）

厚かましいようですが、式の値の正負の記号についてもう一度確認していただけますでしょうか？

よろしくお願いいたします。

回答No.6

info22

#1～#5です。

A#5の補足の質問の回答

>任意定数とは微分方程式の一般解を
> y=C1exp(α1t)+C2exp(α2t)　　・・・(1)
> とした時のC1,C2のことでしょうか？
そうです。
C1,C2とA#4,A#5のD,Eは同じ任意定数です。

> それでしたら
> 初期条件がx(0)=x0(初期位置),dx(0)/dt=v0(初期速度)
> としたときに(1)に代入して
> x0=C1+C2 ・・・(2)
> v0=-α1C1-α2C2　　・・・(3)

>(2),(3)より
> C1=v0/-α1+α2
> C2=v0/α1-α2

> ということでよろしいのでしょうか？
良いですよ。

>その場合，info22さんが求めてくださった，
> f(x)= d*exp((-c+√(c^2-4*a*b))*x/(2*a))+
> e*exp((-c-√(c^2-4*a*b))*x/(2*a))+(-c+x*b)/b^2

> の式のα1,α2（同次方程式の解）を教えていただけるとありがたいです。
d(D)＝C1はα1=(-c+√(c^2-4*a*b))*x/(2*a)に対する任意定数
e(E)＝C2はα2=(-c-√(c^2-4*a*b))*x/(2*a)に対する任意定数
に対応します。

>またX''=(Y-CX'-BX)/Aの式は，
>分類的には定数係数同次線形微分方程式か定数係数非同次線形微分方程式
>のどちらかに分類されると思うのですが，どちらなのでしょうか。

Yの項があるので
定数係数非同次線形微分方程式の方になりますね。

補足

ありがとうございます！
おかげさまでここまでは理解できました。

d(D),e(D)の式のxは#1のxの値に対応しているのですよね？

エクセル(下記参照)でもできるようなのでやってみたいと思います。
http://aoki2.si.gunma-u.ac.jp/Hanasi/StatTalk/solver.html

回答No.5

info22

#1～#4です。
> この場合の未定定数D,Eはどのように扱えばよいのでしょうか？

任意定数C、Dは、2階微分方程式を満たす固定点(x0,y0)における境界条件（初期条件）から決まる定数ですから、

X'(x0,y0)=K1,X(x0,y0)=K2　…(◆)
から決まる定数ですね。

この様な微分方程式の境界条件(◆)がありませんか？
あれば、その条件を適用すればD,Eが決まります。

補足

ご回答ありがとうございます。

任意定数とは微分方程式の一般解を
y=C1exp(α1t)+C2exp(α2t)　　・・・(1)
とした時のC1,C2のことでしょうか？
それでしたら
初期条件がx(0)=x0(初期位置),dx(0)/dt=v0(初期速度)
としたときに(1)に代入して
x0=C1+C2 ・・・(2)
v0=-α1C1-α2C2　　・・・(3)

(2),(3)より
C1=v0/-α1+α2
C2=v0/α1-α2

ということでよろしいのでしょうか？

その場合，info22さんが求めてくださった，
f(x)= d*exp((-c+√(c^2-4*a*b))*x/(2*a))+
e*exp((-c-√(c^2-4*a*b))*x/(2*a))+(-c+x*b)/b^2

の式のα1,α2（同次方程式の解）を教えていただけるとありがたいです。

またX''=(Y-CX'-BX)/Aの式は，
分類的には定数係数同次線形微分方程式か定数係数非同次線形微分方程式
のどちらかに分類されると思うのですが，どちらなのでしょうか。

記号が混同してしまってわかりにくいかもしれませんが
ご返答よろしくお願いします。

回答No.4

info22

#1,#2,#3です。
1週間ほど出かけていて回答が遅れました。
A#3の補足質問について
>X''=(Y-CX'-BX)/Aでフィッティングを行う場合、f(x)=a*x*x*x+b*x*x+c*xとおいてフィッティングを行ってよろしいのでしょうか？
ダメだと思います。
微分方程式を解いた一般解に対してフィッティングして係数を決めます。
一般解の未定定数をD,Eとすると、微分方程式の一般解は次のように求まります。
X=f(Y)= D*exp(1/2*(-C+√(C^2-4*B*A))*Y/A)+
E*exp(-1/2*(C+√(C^2-4*B*A))*Y/A)+(-C+Y*B)/B^2
この関数を改めて
　f(x)= d*exp((-c+√(c^2-4*a*b))*x/(2*a))+
e*exp((-c-√(c^2-4*a*b))*x/(2*a))+(-c+x*b)/b^2
と置いてフィッティングすれば良いでしょう。

この場合、文字定数の項「√(c^2-4*a*b)」が含まれているので改めて
xの係数が単純になるように
(a,b,c,d,e)の文字定数を別のセットに置換えた方が良いかも知れませんね。

補足

ご返答ありがとうございます。

なるほど、微分方程式の一般解に対してフィッティングすればいいのですね。

しかし、この場合の未定定数D,Eはどのように扱えばよいのでしょうか？
Gnuplotでは，C言語などと同じように記号を宣言しなければならないのですが…

何度も何度も申し訳ありません。
おかげさまでもう少しで求まりそうです。
よろしくお願いします。

回答No.3

info22

#1、#2です。

A#2の補足質問について

>最少二乗法による解き方があるみたいなのですが、ご存知でしょうか？

最小二乗法による場合もフィッティング関数を与えてパラメータを求めますので、決定すべき未定係数を含んだフィッテング関数（近似する関数の形）が確定していないと使えません。
今回のように係数によって3通りの関数の形になりますので、データによって、3通りの関数形についてフィってリングを行い、最も適合する関数を人が介在して判断することになるかと思います。

フリーソフト（数式処理および”2D-3Dプロットソフト）の中で、フィッテッィングが特に得意なソフトに「gnuplot」があります。有名なソフトなのでGoogleで調べればダウンロードサイトやフィッティングの日本語の使用例（サンプルプログラム）が沢山ひっかかってきますので、それらを利用されるといいと思います。テキスト形式のx,yのデータファイルを作って予めgnuplotのデレクトリに置いていて、フィッテングコマンドで入力データを読み込み、フィッティング関数を与えておいて、実行すれば、フィッティング関数の係数セットを最小二乗法で求めてくれます。

ダウンロードとインスロール
http://www.geocities.jp/kaname78web/gpdwld.html
Gnuplotでの関数fitting
http://t16web.lanl.gov/Kawano/gnuplot/misc2.html
http://wwwal.kuicr.kyoto-u.ac.jp/www/accelerator/a4/gnuplot/fit.html

補足

返答が遅くなりました、申し訳ありません。

ご回答を参考に、Gnuplotを用いてフィッテングを行ってみました。

X''=(Y-CX'-BX)/Aでフィッティングを行う場合、f(x)=a*x*x*x+b*x*x+c*xとおいてフィッティングを行ってよろしいのでしょうか？

ご返答よろしくお願いします。

回答No.2

info22

#1です。

微分方程式を解くとC1とC2を不定定数として、解くと
C^2>4ABの場合
X=C2*exp[-{C-√(C^2-4AB)}Y/(2A)]+C1*exp[-{C+√(C^2-4AB)}Y/(2A)]
+(YB-C)/(B^2)
C^2<4ABの場合
X=C2*exp[-{C-i√(4AB-C^2)}Y/(2A)]+C1*exp[-{C+i√(4AB-C^2)}Y/(2A)]
=[exp{-CY/(2A)}]〔C1'cos[{√(4AB-C^2)}Y/(2A)]+C2'sin[{√(4AB-C^2)}Y/(2A)]〕
+(YB-C)/(B^2)
C^2=4ABの時
X=(C1"+C2"Y)*exp{-CY/(2A)}+(YB-C)/(B^2)

未知数が5個のA,B,C,C1,C2(C1',C2'; C1",C2")なので
A#1の補足の(X,Y)の５組を与えれば、5変数の連立方程式が得られ
解けるはず。
しかし、A,B,Cによって場合分けがあったり、一次連立方程式ではないので解くのが困難ですね。

微分方程式の解のタイプを絞ったり、変数を減らすなどしないと解くのが難しいかも知れないですね。

補足

ご丁寧なご回答ありがとうございます。

制御に用いたいのですが、場合分けだと困難ですね。
最少二乗法による解き方があるみたいなのですが、ご存知でしょうか？
もし、ご存知でしたらご教授願いたいです。

よろしくお願いします。