分母分子の次数が同じ場合のラプラス逆変換

皆さんよろしくお願いいたします。
次の関数のラプラス逆変換をどのように解いたらよいかわからず困っております。
G(s)=-{(s-a)(s-σ-jω)(s-σ+jω)}/{(s+a)(s+σ-jω)(s+σ+jω)}
ここでjは虚数単位、a,σ,ωは実系数とします。
分母と分子の次数が同じため、次のように部分分数展開しようとすると
分子の次数が合わないため求められません。
G(s)=-(s-a){(s-σ)^2+ω^2}/[(s+a){(s+σ)^2+ω^2}]
=A/(s+a)+{C(s+σ)+Dω}/{(s+σ)^2+ω^2}
この場合どのように部分分数展開し、逆ラプラス変換すればよいのか
ご教示いただければ幸いです。

No.12 ベストアンサー 回答者： stomachman 回答日時： 2008/05/03 01:41

ANo.9のコメントについてです。

> せっかくご教示いただいたのですが、高等すぎてよく分かりません。

ふーむ。γ(n,T,t)のグラフを描いてみてもそう仰るかなあ？ （もしまだでしたら、まずはT=1でn=1～10ぐらいまで描いてみて下さいな。）

　ANo.5では
γ(n,T,t) = (t^(n-1)) exp((1-n)t/T)
を
A(n,T)= ∫γ(n,T,t) dt　　　(積分はt=0～∞）
で規格化したもの
f(n,T,t) = γ(n,T,t)/A(n,T)
を考えた。

γ(n,T,t) = (t^(n-1)) exp((1-n)t/T)
のカーブ下面積は
A(n,T)= ∫γ(n,T,t) dt　　　(積分はt=0～∞）
= ∫(t^(n-1)) exp((1-n)t/T)dt
変数変換
t = xT/(n-1)
で
A(n,T)=((T/(n-1))^n) ∫(x^(n-1)) exp(-x)dx　　(積分はx=0～∞）
となるけど、この積分は「第二種オイラー積分」と呼ばれていて、ガンマ関数を表している。すなわち、
Γ(z) = ∫(x^(z-1)) exp(-x) dx　　(積分はx=0～∞）
そして、
Γ(n) = (n-1)!
である。なので、
A(n,T)=((T/(n-1))^n)(n-1)!
（従って、ANo.5より、fのラプラス変換は
F(n,T,s) = ( ((n-1)/T)/(s+(n-1)/T) )^n
である。すなわちa/(s+a)のn個直列でT=(n-1)/aとなる。）

　かくて
f(n,T,t) =((((n-1)/T)^n)/(n-1)!) (t^(n-1)) exp((1-n)t/T)
と分かった。(n-1)!がヤラシイのでスターリングの公式
n! = √(2πn)　(n^n)exp(-n+θ/(12n))　　（0＜θ＜1)
を（この公式もガンマ関数から導かれる訳ですが）使うと
f(n,T,t) = ((n/√(2π(n-1)))/T) ((t/T)^(n-1)) exp((1-n)(t/T-1)-θ/(12(n-1))
となる（んじゃないかな）。ともあれ、このfは、
・　f(n,T,0)=0
・　t>0のときf(n,T,0)＞0
・　t→∞のとき、f(n,T,0)→0
・　t=Tに唯一の極大を持つ。
・　t=T±T/√(n-1)に変曲点を持つ。
を満たす。

　さて、「任意のtについて、t＞0 ならば、n→∞のときf(n,T,t)→δ(T-t)である」と言いたい。ってことはつまり、
(a) ∫f(n,T,t) dt=1　　　(積分はt=0～∞）
(b)任意のtについて、t＞0かつt≠T ならば　n→∞のときf(n,T,t)→0
を示したい。
　でも(a)はfの定義から自明。あとは(b)を証明すれば良い。要するに極限を計算すれば良いんだし、また「或るt(t＞0, t≠T）について、自然数Nが存在して、n>Nならば f(n,T,t)＞f(n+1,T,t)である」と言えればそれでOK。ですが、やっぱり、めんどくさいばかりで結果が分かってる計算、やる気がしないなあ。なので、自分でやってね。

> どうせ勉強するなら実務で使えるものでないと意味がないと考え、
> 実務想定で勉強中です。

　無駄時間をこうやって構成することって、実務であるんだろうか。
　単に一定時間Tだけの遅れを作るのが目的だとすると、低周波の信号なら、A/D変換してバッファしてD/A変換する、という手で、うんと長いTでも非常に良い精度で出来てしまう。Tが短いのなら表面弾性波デバイスという手もある。そんなんじゃ追いつかないほどの高周波信号を扱うのなら、長いケーブル（や導波管やファイバ）を通してやるだけで無駄時間が発生する。
　となると、むしろ「長いケーブルにアナログ信号を通すことが本来の目的なんだけど、その際に波形が分散して崩れるのをうまく補正したい。つまり、無駄時間要素を通ったのと等価になるようにしたい」というような時にこそ、ご検討なさっている話が活きて来るかも知れませんね。
　
　それにしてもパデ近似のリスト、すげえですね。高次のものの因数分解は数値的にやれば出来るでしょうけど、数式処理ソフトがコマンド一発でやってくれないのかなあ。近頃のソフトは性能が高いが、敷居も値段も高くて一般人にはなかなか手が出ません…

この回答へのお礼

ご回答いただきありがとうございます。
A(n,T)=((T/(n-1))^n)(n-1)!を導いていただいたので
f(n,T,t) =((((n-1)/T)^n)/(n-1)!) (t^(n-1)) exp((1-n)t/T)を描くことが出来ました。
エクセルでT=1,2,5の時について描いてみたところ見事にT=1,2,5にピークが
得られるインパルス応答を再現できていたので驚きました。
（しかしながら、T=1では、n=143以降#NUM!のエラーが出てしまいます。
T=2,5についてもnが大きくなると同じエラーが出ます。エクセルの限界でしょうか？）
証明もありがとうございます。分かりやすかったです。
＞ANo.5より、fのラプラス変換はF(n,T,s) = ( ((n-1)/T)/(s+(n-1)/T) )^nである。
＞すなわちa/(s+a)のn個直列でT=(n-1)/aとなる。）
一次遅れ要素の直列結合で無駄時間要素を近似できることが分かりました。
つまりlim{n→∞}F(n,T,s) =e^(-Ts)ということができるのでしょうか。
もし証明できれば素晴らしいと思いますがいかがでしょうか。
ご意見をいただければ幸いです。
また、ご教示頂いた関数γ(n,T,t) = (t^(n-1)) exp((1-n)t/T)
について詳しく知りたいです。制御工学の書物に載っているのでしょうか。
出典やURLなどありましたら、ご教示いただきたくお願いいたします。

お礼日時：2008/05/10 18:51

No.24 回答者： stomachman 回答日時： 2008/05/21 02:56

　ANo.23のコメントについてです。

> lim{|y|→∞}(1+1/y)^y=e
> と定義できるのはなぜでしょうか。色々とネットで調べてみましたが、

おそらく「定義」という言葉の意味をご存知ないのでしょう。こちらをご参照方→ http://oshiete1.goo.ne.jp/qa43691.html

それはさておき、なんでこんな意味のない計算にこだわっていらっしゃるのか、ますます理解できなくなって来ました。というのは、
(1) F(n,T,t)→exp(-Ts) は（f→δ(t-T)を証明した上でラプラス変換する、というやり方で）既に分かっている。
(2) 「色々とネットで調べてみ」たと仰っている。
　既に答が出てる問題に取り組む事に意味があるとすれば、それは、自力でやるからこそではありませんかね。なのに自力でやらずに「色々とネットで調べて」いらっしゃる。となると、まるで意味のない作業じゃないかなあ。はてさて、一体何をなさりたいのでしょうか。

　ところで、ANo.12でも申し上げたけれども、stomachmanとしては答が出ている極限の計算に興味はないなあ。と言いながらのANo.18以降は、従って、大サービスです。だから、

> よって、F(n,T,s) = (y/(1+y))^(yTs＋1)=（ (1+1/y)^(-yTs) )(1+1/y)^(-1)
> ではないでしょうか。

　はい。いつもの事ながらstomachman、アホ中学生レベルの計算間違いをやらかしました。「ではないでしょうか」も何も、間違ってるものは間違ってるし、正しい計算は誰がやっても正しい。で、そんな当たり前のことを確認してどうなさるお積もりでしょう？

　という訳で、もし、計算を事細かにチェックして貰いたい、(1+1/y)^yの極限の計算も誰かに教えてもらおう、とお考えなら、別途質問をお立てになるのが良いでしょうね。

この回答へのお礼

ご回答頂きありがとうございます。
大元の質問とは別に、数学の知識に乏しい小生のために
色々とご教示頂きありがとうございます。
理解しやすい方がより良いと思い。
「f→δ(t-T)を証明した上でラプラス変換する、というやり方」
よりも「F(n,T,t)→exp(-Ts)」を示すほうが直接的で分かりやすいと
小生なりに感じたので、質問をさせて頂きました。
stomachmanさんにとっては意味の無い計算でも小生にとっては、
今後ご教示いただいたｆを使いこなす上で、しっかり証明を理解
しておきたかったのです。ご面倒をおかけして申し訳有りません。
また、ご興味の無いことにつき合わせてしまい申し訳ありませんでした。

F(n,T,s) = (y/(1+y))^(yTs＋1)=（ (1+1/y)^(-yTs) )(1+1/y)^(-1)
の計算についても、小生が行ったものが正しいのかどうか自信がないため
ご確認いただきたかったのです。

lim{|y|→∞}(1+1/y)^y=eについても、なるべく自力でやろうと
学生時代の教科書を色々調べたのですが、載っていなかったため、
ネットで調べてみることにしました。

>ANo.12でも申し上げたけれども、stomachmanとしては答が出ている
>極限の計算に興味はないなあ。と言いながらのANo.18以降は、
>従って、大サービスです。
大サービス本当にありがとうございます。小生の乏しい知識では出来ないことです。

>(1+1/y)^yの極限の計算も誰かに教えてもらおう、とお考えなら、
>別途質問をお立てになるのが良いでしょうね。
現在のところ、小生の乏しい知識ではyが複素数の場合の(1+1/y)^yの
極限からeを導き出すことが出来ないでおります。
よって新たに質問を立てることにいたします。
ありがとうございました。

お礼日時：2008/05/21 11:48

No.23 回答者： stomachman 回答日時： 2008/05/20 07:45

> ANo.12にてF(n,T,s) = ( ((n-1)/T)/(s+(n-1)/T) )^n

　このFは1次遅れ系をn個直列したものでした。で、それぞれの1次遅れ系
　　((n-1)/T)/(s+(n-1)/T)
の分子と分母をsで割ったのでしょう？ならば

　　((n-1)/(Ts))/(1+(n-1)/(Ts))

になりませんか？って訳で、

　　y = (n-1)/(Ts)

とおいて、

　　F(n,T,s) = (y/(1+y))^(yTs-1) = (1+1/y)( (1+1/y)^(-yTs) )

　あと、こまかいことですけど、ANo.21のコメントで

> lim{y→∞}(1+1/y)^y=eをつかうとx-1→∞において

とお書きですが、「y→∞」だとyは実数だということになる。yは複素数なんで、

　　lim{|y|→∞}(1+1/y)^y=e

でないとおかしい。

この回答へのお礼

ご回答頂きありがとうございます。
長いことつき合わせてしまいすいません。
お時間を頂き大変恐縮で申し訳ありませんが、今しばらくお付き合いください。
>((n-1)/T)/(s+(n-1)/T)の分子と分母をsで割ったのでしょう？
>ならば((n-1)/(Ts))/(1+(n-1)/(Ts))になりませんか？って訳で、
おっしゃるとおりです。そのように変換し(n-1)/(Ts)をeの定義が自分なりに
使いやすいように-xとおきました。
ご教示いただいたように単にy = (n-1)/(Ts)とおいた方が分かりやすいですね。
ただここでy = (n-1)/(Ts)なのでn=yTs+1になると思います。
よって、F(n,T,s) = (y/(1+y))^(yTs＋1)=（ (1+1/y)^(-yTs) )(1+1/y)^(-1)
ではないでしょうか。

また新たな疑問が出てきてしまいました。
> lim{y→∞}(1+1/y)^y=eをつかうとx-1→∞において
>とお書きですが、「y→∞」だとyは実数だということになる。
>yは複素数なんでlim{|y|→∞}(1+1/y)^y=e
ここでeの定義をネットなどで調べたところ、yは実数で定義されています。
複素数のyの場合のeがご教示いただいたように
lim{|y|→∞}(1+1/y)^y=e
と定義できるのはなぜでしょうか。色々とネットで調べてみましたが、複素数
の場合でｅを定義しているところを見つけられませんでした。
数学的にどのように定義して上式を導けるのか、お時間のあるときで結構ですので、
ご存知であればご教示いただきたくお願いいたします。

お礼日時：2008/05/20 22:52

No.22 回答者： stomachman 回答日時： 2008/05/20 00:03

ANo.21のコメントについて。

> ご回答が遅くなり大変申し訳有りません。

どっちが回答者か分からなくなって来ました (^^;
えーと、

> 左辺=lim{n→∞}[( ( (n-1)/(Ts) )/( s+(n-1)/(Ts) )]^n

これ、計算間違いしてませんか？

この回答へのお礼

ご連絡頂きありがとうございます。
ANo.12にてF(n,T,s) = ( ((n-1)/T)/(s+(n-1)/T) )^n
とお示しになられているので上式となります。上式を変換すると以下のようになります。
F(n,T,s) = ( ((n-1)/T)/(s+(n-1)/T) )^n=[( ( (n-1)/(Ts) )/( s+(n-1)/(Ts) )]^n
と変換しました。
まちがっているでしょうか。

お礼日時：2008/05/20 02:46

No.21 回答者： stomachman 回答日時： 2008/05/14 02:41

　ANo.19のコメントについてです。

> なぜF(n,T,s) をラプラス逆変換したf(n,T,t)について
>「t＞0かつt≠Tのとき、lim{n→∞} f(n,T,t) = 0　を示す」ことで
> lim{n→∞}F(n,T,s) =e^(-Ts)を証明したことになるのでしょうか。
> まずは、ここでつまづいてしまいました。

　うーん。どうもこう、なんかすれ違ってますねえ…

　それで証明になっているかどうか、というご質問に関する直接の答は、（既に説明してあるけれども）
　∫f(n,T,t)dt = 1 (積分はt=0～∞）
と分かっていますから、後は
　t≠Tならば lim f(n,T,t) → 0
を示せば、
　f(n,T,t)→δ(t-T)
が言えたことになる。（ラプラス変換の方は、f→δ(t-T)なんだから、F→e^(-Ts)も当たり前。）

　もちろん、F(n,T,t)→exp(-Ts)の証明をやっても良い。公式
　x→∞のとき　 (1-1/x)^(-x)→ e
を使えば出来るんじゃないかな。

　しかしですね、f→δ(t-T) やF→exp(-Ts)を証明することは主要な目的じゃありません。（だって自明だもん。）んな事よりも、（実践を考えるのであれば）nが大きくなるに連れてどんな風に収束して行くかを把握することが重要でしょう。で、それを見るにはラプラス変換Fより、実体であるfの振る舞いを調べる方が適しているだろうと。（実際、パデ近似の方も、大変な苦労をして逆ラプラス変換をやってるのは何のためかというと、「構成した伝達関数の実体はどうなっていて、どんな誤差が残るかを見たいがため」でしょ？）

この回答への補足

ご回答が遅くなり大変申し訳有りません。
F(n,T,t)→exp(-Ts)の証明をやってみましたので、合っているかご教示いただければ幸いです。
lim{n→∞}[( ((n-1)/T)/(s+(n-1)/T) )]^n=e^(-Ts)

左辺=lim{n→∞}[( ( (n-1)/(Ts) )/( s+(n-1)/(Ts) )]^n
(n-1)/(Ts)=-xとおくと
[( ( (n-1)/(Ts) )/( s+(n-1)/(Ts) )]^n=( -x/(1-x) )^( (-Ts)x+1 )
={( -x/(1-x) )^(-Ts)x}( -x/(1-x) )
={( 1+1/(x-1) )^(-Ts)(x-1)}{( 1+1/(x-1) )^(-Ts)}( 1-1/(1-x) )
ここで
lim{y→∞}(1+1/y)^y=eをつかうとx-1→∞において
{( 1+1/(x-1) )^(-Ts)}→1、( 1-1/(1-x) )→1、
{( 1+1/(x-1) )^(-Ts)(x-1)}→e^(-Ts)
よって
lim{n→∞}[( ((n-1)/T)/(s+(n-1)/T) )]^n=e^(-Ts)
が示せた。
いかがでしょうか、間違っている点などご指摘を頂ければ幸いです。

補足日時：2008/05/19 18:06

この回答へのお礼

ご回答頂きありがとうございます。
小生、相変わらず体調不良で寝込んでいます。
返信遅くなり申し訳有りません。
　f(n,T,t)→δ(t-T)
を証明するために、δ関数の性質である。
∫f(n,T,t)dt = 1 (積分はt=0～∞）とt≠Tならば lim f(n,T,t) → 0
を証明したのですね。これで理解できました。ありがとうございます。
小生の理解不足でお手数をおかけして大変申し訳有りません。
＞もちろん、F(n,T,t)→exp(-Ts)の証明をやっても良い。公式
＞　x→∞のとき　 (1-1/x)^(-x)→ e
＞を使えば出来るんじゃないかな。
上の方も納得できましたが、こちらの方が直接的で分かりやすい感じました。
証明にトライしてみたいと思います。今しばらくお時間をください。
＞f→δ(t-T) やF→exp(-Ts)を証明することは主要な目的じゃありません。（だって自明だもん。）
すいません。自明であることが、小生の理解不足、経験不足のため理解できてません。
＞んな事よりも、（実践を考えるのであれば）nが大きくなるに連れて
＞どんな風に収束して行くかを把握することが重要でしょう。
おっしゃるととおりです。しかし将来自分も教える立場になったとき、
stomachmanさんのように、懇切丁寧で分かりやすい説明ができるように
なりたいと思います。
やはり自分だけで納得しているだけでは、社会のために活かせないと思います。
ここの「教えて!goo」もそんな趣旨があると思いました。
いずれ世のため人のためになるよう、今しっかり勉強していきたいと思います。

お礼日時：2008/05/15 13:44

No.20 回答者： stomachman 回答日時： 2008/05/12 15:33

ANo.19のつづきです。

んもお、宙ぶらりんは気持ち悪いから
　n→∞のときf(n,T,t)→0
の証明の流れを書いちゃいます。

　ln(f(n,T,t))においてθ=0にしたものをgとする。（これが-∞に発散してくれれば、f(n,T,t)が0に収束する訳です。）すなわち
　g(n) = ln((n/√(2π(n-1)))/T)+ ln((t/T)^(n-1)) + (1-n)(t/T-1)
　= ln((n/√(2π(n-1)))/T)+ (n-1)ln(t/T) + (1-n)(t/T-1)
　=ln(n)-(1/2)ln(n-1)+lnπ-lnT+(n-1)(ln(t/T)+1-t/T)
ここにt≠T, t＞0, T＞0。

　さて、nはn＞1を満たす実数だと思うことにする。g(n)が任意のn(n＞1)において連続である、ということは自明。 gのnによる1階微分を考える。
　g'(n) = 1/n -1/(2(n-1))+(ln(t/T)+1-t/T)
　　= (n-2)/(2n(n-1))+(ln(t/T)+1-t/T)

(1)もし全てのn（n＞1)についてg'(n)＜0であれば、
　n→∞のとき、g(n)→-∞
である。終わり。

(2)では、もしg'(n)≧0となるn(n＞1)が存在するならばどうなるか。
　h(n)=(n-2)/(2n(n-1))
と定数
　C=ln(t/T)+1-t/T
を使って、
　g'(n)=h(n)+C
と書ける。
　x=1-t/T
とおけば、t＞0とT＞0よりt/T＞0なので、
　x＜1
である。そこでテイラー展開の公式
　ln(1-x)= -Σ((x^k)/k) （Σはk=1～∞、x＜1)
を使うと
　ln(1-x)+x ＜ 0
すなわち
　C ＜ 0
であることが分かる。
　一方、hの極値を求めると、
　h'(n) = -2(n^2+4n-2)/(2n(n-1))^2
より、
　h'(p)=0
を満たすpは p＞1 の範囲では
　p = 2+√2 ≒ 3.4
だけであって、hはここで極大である。従って、
　n＞pのとき、h'(n)＜0
である。
　g'(n)≧0となるn(n＞1)が存在すると仮定したのだから、
　g'(x)=h(x)+C=0
を満たすxが存在する。（これはただの二次方程式であるから、実際に解いてNをCの関数として表す事もできる。）すると、そのようなxは高々２個存在し、そのうち最大のものはx≧pでなくてはならない。それをNと書くことにする。すると、
　n＞N ならば g'(n) ＜ 0
が成り立つのは自明。だから、
　n→∞のとき、g(n)→-∞
終わり。

No.19 回答者： stomachman 回答日時： 2008/05/11 17:51

ANo.18のコメントについてです。

> 対数を取って整理したもの？最後にexp()で変換？

　ANo.12の

>> f(n,T,t) = ((n/√(2π(n-1)))/T) ((t/T)^(n-1)) exp((1-n)(t/T-1)-θ/(12(n-1))

を使います。nが大きいときのだいたいの格好を知りたい訳ですから、θは0にしちゃえばいいですね。

　f(n,T,t) = ((n/√(2π(n-1)))/T) ((t/T)^(n-1)) exp((1-n)(t/T-1))

の両辺の自然対数 ln( ) をとると、

　ln(f(n,T,t)) = ln((n/√(2π(n-1)))/T)+ ln((t/T)^(n-1)) + (1-n)(t/T-1)
　= ln((n/√(2π(n-1)))/T)+ (n-1)ln(t/T) + (1-n)(t/T-1)

となる。Excelで右辺を各tについて計算した数値をX(t)とすると

　f(n,T,t) = exp(X(t))

です。（片対数グラフに慣れていらっしゃるのなら、X(t)のままプロットしたって構わん訳ですが。）


> ∀や∃はどこかで診た覚えがあるので、まずはそこからチャレンジしてみます。

　そういえば見かけます∃ね…(^∀^)。って冗談はさておき、制御工学とはもはや全然関係ない、ただの極限の計算です。論理式で書きましたけど、却って分かりにくかったですか。ならば、もっとなじみがありそうな書き方をすれば

「t＞0かつt≠Tのとき、
lim f(n,T,t) = 0
n→∞
を示せ」

です。t,Tを固定してfをnだけの関数だと思えばいいのです。t＞Tとt＜Tとに分けて考えた方が良いかもしれません。また、nを整数に限定するのではなく、実数だと思って、fをnの連続関数だと考えると良いかも。ところで、f(n,T,t)>0は自明なので、

lim ln(f(n,T,t)) = -∞
n→∞

を示しても同じです。さらに、θ＞0, n＞1であるから、

0 ＜ exp((1-n)(t/T-1)-θ/(12(n-1)) ＜ exp((1-n)(t/T-1))

つまり冒頭に挙げた式において、θ=0を代入したもので証明すれば十分。これでだいぶ簡単になったんじゃないかな。
　この関数は、nが小さいうちはnの増加につれて増加することもある。つまり、nの関数として見た時、極大が高々１個あって、その先では単調減少するんです。ですから、極大のありかを実際に計算してみて、その先を考える、というのもひとつの手でしょう。


　それで思い出したのだけど、ANo.18でガンマ型曲線の話をしましたが、その一般形は

K(t^α)exp(-t/β)

であって、K, α, βは実数。αが整数という縛りはこの場合、ないんです。α≦0やβ≦0の場合すら使われることがあります。

この回答へのお礼

体調を崩して寝込んでいるmathstudyです。
レス遅くて申し訳有りません。
ご回答いただきありがとうございます。
f(n,T,t) =((((n-1)/T)^n)/(n-1)!) (t^(n-1)) exp((1-n)t/T)
についてのエクセルで描かせるためのテクニックとして
式の対数を取って整理したものを計算し、最後にexp( )で変換する方法
については、懇切丁寧なご回答を頂いたおかげで理解できました。
本当にありがとうございます。

∀ε(ε>0 ⇒ ∃N∀n(n＞N ⇒ |f(n,T,t)|<ε))
の論理式の意味についても、分かりやすく示していただきありがとうございます。
また証明についても、懇切丁寧にお示し頂きありがとうございます。
こちらの方については、まだ体調不良で取り掛かれておりません。
小生の都合で申し訳有りませんが、今しばらくお時間を下さい。

ところで、小生の理解不足で申し訳有りませんが、
lim{n→∞}F(n,T,s) =lim{n→∞}[( ((n-1)/T)/(s+(n-1)/T) )]^n=e^(-Ts)
を証明するのに、なぜF(n,T,s) をラプラス逆変換したf(n,T,t)について
「t＞0かつt≠Tのとき、lim{n→∞} f(n,T,t) = 0　を示す」ことで
lim{n→∞}F(n,T,s) =e^(-Ts)を証明したことになるのでしょうか。
まずは、ここでつまづいてしまいました。
小生の理解不足（数学音痴）で申し訳ありませんが、ご教示いただければ幸いです。

お礼日時：2008/05/13 11:58

No.18 回答者： stomachman 回答日時： 2008/05/10 21:32

ANo.12のコメントについてです。

> エクセルの限界でしょうか？

　n=143が限界ってのは、おそらく式をそのまま計算なさったのでしょうね。式の対数を取って整理したものを計算し、最後にexp( )で変換してやれば、かなり大きいnまで行けると思います。

> つまりlim{n→∞}F(n,T,s) =e^(-Ts)ということができるのでしょうか。

　証明はご自分でできる程度のものだと思いますよ。

>> (b)任意のtについて、t＞0かつt≠T ならば　n→∞のときf(n,T,t)→0

　つまり、t＞0かつt≠T のとき、

∀ε(ε>0 ⇒ ∃N∀n(n＞N ⇒ |f(n,T,t)|<ε))

を示せば良いのです。

> また、ご教示頂いた関数γ(n,T,t) = (t^(n-1)) exp((1-n)t/T)
> について詳しく知りたいです。制御工学の書物に載っているのでしょうか。

　ANo.5の最後に書いたラプラス変換（ただし分子が１のもの）なら載っているのではないでしょうか。この関数は数学の用語で言えば「不完全ガンマ関数の1階導関数」に過ぎないのですが、いろんな現象に当て嵌まることが多いので、データをパラメータで簡潔に表現して整理するのに使われます。
　また、系がこの関数で表されることが理論的に自明でない応用の例として、血液が肺の血管網（無数に分岐した血管）を通過する時間の分布がこの関数で近似できる、という話があり、その分野ではgamma-shaped curve（ガンマ型曲線）と呼ばれることもあります。
　stomachman的には、「微小な一次遅れ系が無数に直列したものは、無駄時間系になる」ということに気付けたのが大収穫でした。たとえば、弾性のあるパイプにインパルスで圧力を送ると、波形が崩れずに向こう側に到達する、ってことですよね。

この回答へのお礼

ご回答いただきありがとうございます。
＞n=143が限界ってのは、おそらく式をそのまま計算なさったのでしょうね。
そのとおりです。短絡思考で申し訳有りません。
＞式の対数を取って整理したものを計算し、最後にexp( )で変換してやれば、かなり大きいnまで行けると思います。
対数を取って整理したもの？最後にexp()で変換？
とりあえずエクセルを色々いじくってチャレンジしてみます。
よく原理が理解できません。小生の理解不足で申し訳有りません。
もしお時間がおありでしたら具体的にどのようにすればよいかご教示いただければ幸いです。
＞lim{n→∞}F(n,T,s) =e^(-Ts)証明はご自分でできる程度のものだと思いますよ。
すいません。小生統計学を学んでいたことが有り、Γ関数やスターリングの公式
についてはある程度予備知識があったので、すんなり理解できましたが、
制御工学については全くの初学者で、数学の嵐に少々参り気味です。
＞つまり、t＞0かつt≠T のとき、∀ε(ε>0 ⇒ ∃N∀n(n＞N ⇒ |f(n,T,t)|<ε))を示せば良いのです。
申し訳有りません。数学の記号からネットで調べて勉強してみます。
∀や∃はどこかで診た覚えがあるので、まずはそこからチャレンジしてみます。
ほんとうにお時間のあるときで結構ですので、証明の節目節目だけでも
よろしかったら、ご教示いただければ幸いです。
小生の知識に乏しいため、勉強してからお礼を書き込みするまで時間がかかり大変申し訳有りません。
ご回答を、直ぐにいただけているのに、本当に申し訳有りませんが、
しばらくお時間を頂きたくお願いいたします。

お礼日時：2008/05/11 15:40

No.17 回答者： info22 回答日時： 2008/05/08 04:11

#2,#6,#10,#11,#13,#14,#15,#16です。

A#16の補足質問の解答です。

無料ソフトには利用者の要求を100%満たすものはないということを念頭において下さい。その上で上手く使いこなす事です。
数式処理ソフトで一番優先されるべきことは、計算精度を優先する事。次に使用者が使いやすいように余分な入力や操作を減らしできるだけ自動計算してくれること。色々な計算に対応できる事。必要なだけ精度が取れること。
といったことですね。
他の不便さは、ソフトの計算機能を有効に使ってカバーすること。これ以上要求しても、自動的に計算して計算精度が落ちていた、では困ります。
不満足であれば、数式処理ソフトを使って係数等を計算させたり、内蔵の関数がある場合は利用して計算させる。といった方法で我慢する事ですね。無償提供の相手に完全な機能を要求できません。
といっても、有料の数式処理ソフトでも完全ではないですね。人が手助けをしながら、結果を導いたり、不完全な所は、人の能力で補って完全な結果を得る事が大切です。数式処理ソフトに要求すべき機能は、上述の計算精度を必要なだけ計算できて、できるだけ人の手を煩わせない自動化計算ができることです。そして無償ソフト開発している人たちによって、毎年のようにバージョンアップされ、新たな計算機能と不十分な機能の改善がなされてソフトであることです。

これらのことを踏まえた上で
補足質問 1)について

> (%i6) f3:float(ilt(f2,s,t));を実行したら、
> 指数がすべて%eで表示されずに実数表示2.718281828459045になってしまう。
f3:float(ilt(f2,s,t));
と
f3:ilt(f2,s,t);
どちらがよいか、より良いを使って下さい。

係数の有効桁数は利用者の方で減らして下さい。
係数は数式処理ソフトで個別に計算して必要な精度と式にまとめる。
部分計算ではsimplify関数やfactor関数を使うと良い。

補足質問2)について
> 有限の桁数にならず以下のように有効数字のない小数点の指数表示になる。すべて%eで表示されずに実数表示2.718281828459045になってしまう。
これはfloat関数の機能がe=2.71828...とexp(...)関数で上手く使い分ける機能がMaximaではまだ備わっていないから仕方がないですね。exp(...)のexpとeをMaximaのfloat関数が同一として数値化してしまっているのでしょう。人が判断すればすむ事です。

できるだけ計算精度を落さないようにして、整数係数にする為、10^n　の桁移動調整をしているのだと思います。勝手に桁落ち計算されるよりは余程ましです。
汎用プログラム言語（FortranやCによる数値計算）による計算では見かけ上の有効桁数だけ確保するため、実際の精度が１～２桁になっていても何も表示せず見かけ上の決められた桁数で計算結果を出します。数式処理ソフトは計算精度を優先します。

係数が多桁になった場合は使う人の判断で有効桁数をカットしましょう。

> これ自体はmaximaの正確さを物語っているのでしょうが、ご回答頂いた
様にすっきりした形にならず扱いに困っています。

必要な解答を得るには、人の判断による式の整理やソフトによる係数の計算を行います。人が望むすっきりした式にならない場合は、人の手で式を整理したりします。係数等の計算は数式処理ソフトを使って行うのが当然です。

という事で、数式処理ソフトをうまく使いこなす事です。不完全な所は人が手助けして望む結果を得るようにして下さい。

他の有料の数式処理ソフトのMapleで同じ計算をすると、
solve(15120+8400*s+2100*s^2+300*s^3+25*s^4+s^5 = 0, {s})；
では
{s = RootOf(15120+8400*_Z+2100*_Z^2+300*_Z^3+25*_Z^4+_Z^5, index = 1)}, {s = RootOf(15120+8400*_Z+2100*_Z^2+300*_Z^3+25*_Z^4+_Z^5, index = 2)}, {s = RootOf(15120+8400*_Z+2100*_Z^2+300*_Z^3+25*_Z^4+_Z^5, index = 3)}, {s = RootOf(15120+8400*_Z+2100*_Z^2+300*_Z^3+25*_Z^4+_Z^5, index = 4)}, {s = RootOf(15120+8400*_Z+2100*_Z^2+300*_Z^3+25*_Z^4+_Z^5, index = 5)}
という解しか出ません。
これを
evalf(solve(15120+8400*s+2100*s^2+300*s^3+25*s^4+s^5 = 0, {s}))；
とすることで
{s = -3.655694325+6.543736899*I}, {s = -5.700953299+3.210265600*I}, {s = -6.286704752}, {s = -5.700953299-3.210265600*I}, {s = -3.655694325-6.543736899*I}
という結果に変わります。

ラプラス逆変換も
> inttrans[invlaplace](f4, s, t);
では、不要な10^n項が表示されますので、
> simplify((inttrans[invlaplace])(f4, s, t))
とすることで
273.3432917*exp(-6.286704728*t)
+31.65360390*exp(-3.655694321*t)*cos(6.543736899*t)
+48.25129172*exp(-3.655694321*t)*sin(6.543736899*t)
-299.9968954*exp(-5.700953315*t)*cos(3.210265600*t)
-136.0845635*exp(-5.700953315*t)*sin(3.210265600*t)
とできますが、
exp(-3.655694321*t)やexp(-5.700953315*t)の項がまだ括りだせていませんね。(Maximaも同じ。)
Maximaよりは整理された式になっています。
望む式の形にするには人の手を加えてやらないといけませんね。
f(t)=
273.3432917*exp(-6.286704728*t)
+(31.65360390*cos(6.543736899*t)+48.25129172*sin(6.543736899*t))
*exp(-3.655694321*t)
-(299.9968954*cos(3.210265600*t)+136.0845635*sin(3.210265600*t))
*exp(-5.700953315*t)

数式処理ソフトに過大な期待をし過ぎるのは無理があります。上手く使いこなせば便利で役立つソフトです。何よりも良い特長は精度の良い計算ができることです。また絶えず進化して改善されていくソフトであるということです。

No.16 回答者： info22 回答日時： 2008/05/07 20:17

#2,#6,#10,#11,#13,#14,#15です。

A#15の補足質問の回答
>分母の解を求め分母を一次式(s-a)、
>共役複素数の2次式{(s-σ)^2+ω^2}=s^2-2σs+(ω^2+σ^2)の形に
>因数分解し、分子に4次式を代入してparfrac関数で
>部分分数展開を試みましたが、うまくいきません。

A#15の回答の内容をそのままMaximaでやっていただけば出来るはずですが
全然理解頂けていないようですので
Maximaでの全入力、全出力を全部貼り付けておきます。
ただし、出力の余分な出力や必要以上については
桁数表示を減らしたり、自然対数の底を「%e」で置換して
ここでの回答文書を短くしました。

以下Maximaの全入出力（pade近似:(4次)/(5次),L=1の場合)です。

(%i1) f:taylor(exp(-s),s,0,9);f1:pade(f,4,5);
(%o1) 1-s+s^2/2-s^3/6+s^4/24-s^5/120+s^6/720-s^7/5040+s^8/40320
-s^9/362880+...
(%o2) [(5*s^4-120*s^3+1260*s^2-6720*s+15120)
/(s^5+25*s^4+300*s^3+2100*s^2+8400*s+15120)]
(%i3) allroots(s^5+25*s^4+300*s^3+2100*s^2+8400*s+15120);
(%o3) [s=3.210265600308549*%i-5.70095329867182,
s=-3.210265600308549*%i-5.70095329867182,
s=6.543736899360078*%i-3.655694325463562,
s=-6.543736899360078*%i-3.655694325463562,s=-6.286704751729237]
(%i4) expand((s+5.70095329867182)^2+3.210265600308549^2)
*expand((s+3.655694325463562)^2+6.543736899360078^2)
*(s+6.286704751729237);
(%o4) (s+6.286704751729237)
*(s^2+7.311388650927124*s+56.18459360927314)
*(s^2+11.40190659734364*s+42.80667373816151)
(%i5) f2:float(partfrac((5*s^4-120*s^3+1260*s^2-6720*s+15120)
/((s+6.286704751729237)
*(s^2+7.311388650927124*s+56.18459360927314)
*(s^2+11.40190659734364*s+42.80667373816151)),s));
(%o5) -(2.027521869620528*10^8*(2.5625850880652117*s+18.340917962065043))
/(1731917.0*s^2+1.9747156*10^+7*s+7.4137607*10^+7)
+(4.4004918436100836*10^7*(1.0298702344625587*s+14.037814970170099))
/(1431728.0*s^2+1.046792*10^+7*s+8.0441055*10^+7)
+670237.7538028954/(2452.0*s+15415.0)
(%i6) f3:float(ilt(f2,s,t));
(%o6) (6.9845669009756039*10^-7
*(6.9082727672469079*10^+7*sin(6.543736825399252*t)
+4.5319355667292945*10^+7*cos(6.543736825399252*t)))
*%e^(-(15215*t)/4162)
+(5.7739487515856702*10^-7
*(-2.3568717510409307*10^+8*sin(3.210265626289453*t)
-5.1956973088156629*10^+8*cos(3.210265626289453*t)))
*%e^(-(11362*t)/1993)
+273.3432911908646*%e^(-(15415*t)/2452)

この回答へのお礼

ご回答頂きありがとうございます。
全入出力を記載頂きご面倒をおかけして申し訳有りません。
実行してみましたが、以下の問題点が有り困っています。
ご教示いただければ幸いです。

1)(%i6) f3:float(ilt(f2,s,t));を実行したら、
指数がすべて%eで表示されずに実数表示2.718281828459045になってしまう。

2)(%i5) f2:float(partfrac((5*s^4-120*s^3+1260*s^2-6720*s+15120)/((s+6.286704751729237)*(s^2+7.311388650927124*s+56.18459360927314)*(s^2+7.311388650927124*s+56.18459360927314)*(s^2+11.40190659734364*s+42.80667373816151)),s));
を実行したら、(%o5)に示されたように有限の桁数にならず
以下のように有効数字のない小数点の指数表示になる。
-(2.1706258313418866*10^-69*(6.7207667854674306*10^+75*s+4.0092825057882691*10^+76))/(1731917.0*s^2+1.9747156*10^+7*s+7.4137607*10^+7)+(1.0562360981279383*10^-80*(3.9689706743263774*10^+86*s+1.377548962389403*10^+87))/(1431728.0*s^2+1.046792*10^+7*s+8.0441055*10^+7)+(4.4004918436100836*10^-40*(1.4744940510466105*10^+53*s+2.0098332751611698*10^+54))/(1431728.0*s^2+1.046792*10^+7*s+8.0441055*10^+7)^2+13474.09082919691/(2452.0*s+15415.0)
これ自体はmaximaの正確さを物語っているのでしょうが、ご回答頂いた
様にすっきりした形にならず扱いに困っています。

お礼日時：2008/05/07 22:28