こんにちはゲストさん。会員登録(無料)して質問・回答してみよう!

解決済みの質問

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

皆さんよろしくお願いいたします。
次の関数のラプラス逆変換をどのように解いたらよいかわからず困っております。
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}
この場合どのように部分分数展開し、逆ラプラス変換すればよいのか
ご教示いただければ幸いです。

投稿日時 - 2008-04-28 18:08:52

QNo.3982668

困ってます

質問者が選んだベストアンサー

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

投稿日時 - 2008-05-03 01:41:23

お礼

ご回答いただきありがとうございます。
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:27

このQ&Aは役に立ちましたか?

2人が「このQ&Aが役に立った」と投票しています

回答(24)

ANo.24

 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の極限の計算も誰かに教えてもらおう、とお考えなら、別途質問をお立てになるのが良いでしょうね。

投稿日時 - 2008-05-21 02:56:49

お礼

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

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:01

ANo.23

> 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

でないとおかしい。

投稿日時 - 2008-05-20 07:45:39

お礼

ご回答頂きありがとうございます。
長いことつき合わせてしまいすいません。
お時間を頂き大変恐縮で申し訳ありませんが、今しばらくお付き合いください。
>((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
と定義できるのはなぜでしょうか。色々とネットで調べてみましたが、複素数
の場合でeを定義しているところを見つけられませんでした。
数学的にどのように定義して上式を導けるのか、お時間のあるときで結構ですので、
ご存知であればご教示いただきたくお願いいたします。

投稿日時 - 2008-05-20 22:52:59

ANo.22

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

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

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

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

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

投稿日時 - 2008-05-20 00:03:55

お礼

ご連絡頂きありがとうございます。
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:58

ANo.21

 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の振る舞いを調べる方が適しているだろうと。(実際、パデ近似の方も、大変な苦労をして逆ラプラス変換をやってるのは何のためかというと、「構成した伝達関数の実体はどうなっていて、どんな誤差が残るかを見たいがため」でしょ?)

投稿日時 - 2008-05-14 02:41:34

補足

ご回答が遅くなり大変申し訳有りません。
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:46

お礼

ご回答頂きありがとうございます。
小生、相変わらず体調不良で寝込んでいます。
返信遅くなり申し訳有りません。
 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:03

ANo.20

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は高々2個存在し、そのうち最大のものはx≧pでなくてはならない。それをNと書くことにする。すると、
 n>N ならば g'(n) < 0
が成り立つのは自明。だから、
 n→∞のとき、g(n)→-∞
終わり。

投稿日時 - 2008-05-12 15:33:59

ANo.19

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の関数として見た時、極大が高々1個あって、その先では単調減少するんです。ですから、極大のありかを実際に計算してみて、その先を考える、というのもひとつの手でしょう。


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

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

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

投稿日時 - 2008-05-11 17:51:36

お礼

体調を崩して寝込んでいる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:19

ANo.18

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

投稿日時 - 2008-05-10 21:32:57

お礼

ご回答いただきありがとうございます。
>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:21

ANo.17

#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による数値計算)による計算では見かけ上の有効桁数だけ確保するため、実際の精度が1~2桁になっていても何も表示せず見かけ上の決められた桁数で計算結果を出します。数式処理ソフトは計算精度を優先します。

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

> これ自体は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)

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

投稿日時 - 2008-05-08 04:11:57

ANo.16

#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)

投稿日時 - 2008-05-07 20:17:38

お礼

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

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:26

ANo.15

#2,#6,#10,#11,#13,#14です。
A#14の補足質問の回答

>factor関数による因数分解や、partfrac関数による部分分数展開を試みましたが、
うまくいきません。
factor関数は有理係数の範囲での因数分解です。
parfrac関数は有理係数の多項式因数での部分分数展開ですから、factor関数で因数分解できない多項式が分母にくる場合は因数分解できる(できている)分母因数について部分分数展開できると思います。分母を一次式(s-a)、共役複素数の2次式{(s-σ)^2+ω^2}=s^2-2σs+(ω^2+σ^2)の形に因数分解しておいてやれば良いかと思います。

Maximaの場合
分母=0が複素解を持つ場合、数値計算で全ての解を求めるには
algsys([f(s)=0],[s]);
または
allroots(f(s)=0);
のどちらの関数でも求まります。

なお、Mapleでは solve関数で複素解をもつ場合も解けます。

投稿日時 - 2008-05-07 10:33:31

お礼

ご解答頂きありがとうございます。
ご教示頂いた関数で、分母の解を求め分母を一次式(s-a)、
共役複素数の2次式{(s-σ)^2+ω^2}=s^2-2σs+(ω^2+σ^2)の
形に因数分解し、分子に4次式を代入してparfrac関数で
部分分数展開を試みましたが、うまくいきません。
無理やり上記の状態でilt関数でラプラス逆変換しましたが、うまくいきません。
本当に質問から外れたものになって申し訳ありませんが、
maximaでのやり方をご存知でしたら、ご教示いただきたくお願いいたします。

投稿日時 - 2008-05-07 17:43:19

ANo.14

#2,#6,#10,#11,#13です。
> 5次以上の高次の解を求めるには、何かテクニックがあるのでしょうか。
L=1,(4次)/(5次)のパデ近似についてもやって見ました。
分母が5次でも数値計算で因数分解でき
exp(-s)~{1-(4/9)*s+(1/12)*s^2-(1/126)*s^3+(1/3024)*s^4}
/{1+(5/9)*s+(5/36)*s^2+(5/252)*s^3+(5/3024)*s^4+(1/15120)*s^5}
分母=(1/15120)*(s+6.286704752)
*(s^2+7.311388650*s+56.18459360)*(s^2+11.40190660*s+42.80667374)
Laplace逆変換して
f(t)=273.3432899*exp(-6.286704752*t)
+31.65360373*exp(-3.655694325*t)*cos(6.543736899*t)
+48.25129152*exp(-3.655694325*t)*sin(6.543736899*t)
-299.9968937*exp(-5.700953300*t)*cos(3.210265598*t)
-136.0845591*exp(-5.700953300*t)*sin(3.210265598*t)

となります。
プロットするとt=1(L=1)にピークが出ますが、t<1でまだかなり振動しますね。

同じくL=1の場合の(5次)/(6次),(6次)/(7次),
(7次)/(8次),(8次)/(9次)のパデ近似についても分母の因数分解をやってみましたが、数値計算で容易に因数分解ができました。

一般の5次以上の方程式は解けないですが、
係数が有理数のn次多項式=0の方程式は数値計算的に解けるようですね。(一般解の公式は難しいかも知れませんね。)

分母が(s-a)形の因数と{(s-σ)^2+ω^2}形の項の因数の積に因数分解できることから、逆ラプラス変換も数式処理ソフトで容易に求まります。

パデ近似関数の、分母の次数より低い場合について、色々な次数のパデ近似関数についてもδ(t-L)どの位再現できるか、やってみると良いかと思います。

数式処理ソフトのMathmatica、Mapleは多くの大学や高専にライセンス導入されていて大学生や院生なら無料で使えると思います。Mathematicaはスチューデント版が出ていますので学生なら3万円前後で個人で入手できると思います。
またMaximaはフリーソフトですので誰でもネットから入手して使えます。無料ソフトですので非常に長い式などは扱えない制約がありますが、それでも他の有料ソフト並みのことができます。Microsoft Excelよりずっと使い物になりますので使ってみて下さい。

手計算では時間的な制約が多いので、ぜひ数式処理ソフトを使用してみてください。

投稿日時 - 2008-05-05 18:11:21

お礼

ご回答頂きありがとうございます。
MathmaticaやMapleは高額すぎて、ちょっと個人では入手困難ですが、
Maximaを入手してみました。その能力にちょっと驚きです。
>手計算では時間的な制約が多いので、ぜひ数式処理ソフトを使用してみてください。
全くそのとおりです。有用なフリーソフトをご紹介頂き大変ありがとうございます。
パデ近似もご教示頂いたやり方で簡単に求めることが出来て感激です。
そのほかの方法については、まだマニュアルと苦戦中です。
ご回答頂いた式を導くことさえできずにいます。
本来の質問とはかけ離れ申し訳ありませんが、もしMaximaのコマンドでの
やり方をご存知でしたらご教示いただきたくお願いいたします。
factor関数による因数分解や、partfrac関数による部分分数展開を試みましたが、
うまくいきません。

投稿日時 - 2008-05-05 23:16:36

ANo.13

#2,#6,#10,#11です。

> 小生3次方程式までしか勉強しておりません。
> 5次以上の高次の解を求めるには、何かテクニックがあるのでしょうか。

計算が大変そうな様子でしたので補足しておきます。

パデが(3次)/(4次)の場合が以下のようにできます。
分母が(5次)以上の場合も分母の1次の因数は
数値計算で簡単に求められますのでそれらを全部分離すればいいですね。
そうすると残りの因数は2次因数の積が残ります。
残りの因数についても共役極の2次式因数を1ずつ数値計算的に抜き出して行けば良いかと思います。(これは回路合成のFoster展開を求めることに相当します)

A#10で挙げた(3次)/(4次)のパデの L=1について逆変換を求めてみました。
数値計算でパデ関数の極の1つを求め、その共役複素数を使って分母の2次因数を求め、分母を強引に因数分解してみると
(s^2+6.425613794*s+33.10449180)*(s^2+9.574386206*s+25.37420013)
となり、ラプラス逆変換すると
f(t)=175.5744076*exp(-3.212806897*t)
*{0.1287378970*cos(4.773087433*t)-0.1420671273*sin(4.773087433*t)}
+162.8013340*exp(-4.787193103*t)
*{-0.1634082433*cos(1.567476419*t)+0.7379759268*sin(1.567476419*t)}
f(t)をプロットすると
t=L=1で正のピークが出ていますがそれ以下のtで振動がありますね、
t>2.5では殆どf(t)はゼロに張り付きます。

分母が(4次)のパデは上のようにできますが
分母が(5次)のパデは実数の極を持ちますので、その極を数値計算で求めて分離すれば、分母は(4次)になりますので上記の方法で(5次)のパデについても逆変換f(t)を求めるのは難しいことでは無いと思います。

投稿日時 - 2008-05-03 11:39:47

ANo.11

#2,#6,#10です。
#10の補足質問の回答です。

> どうやったらこれだけ多くの計算を短時間に出来たのでしょうか?

pade近似関数は数式処理ソフトのMaple,Mathmatica,Maximaなどを使えば直ぐ作成します。
例えば、Mapleでは次のようにプログラムを組めばいいですね。
以下は分子が3次、分母が4次のパデ近似を求める例です。
> with(numapprox):
pade(exp(-L*s), s, [3,4]);
結果は
(-4*L^3*s^3+60*L^2*s^2-360*L*s+840)
/(L^4*s^4+16*L^3*s^3+120*L^2*s^2+480*L*s+840)   (1)
となります。

同じことをMaximaですると次のようにプログラムをします。
(%i1) taylor(exp(-L*s),s,0,7);pade(%,3,4);
(%o1) 1-L*s+(L^2*s^2)/2-(L^3*s^3)/6+(L^4*s^4)/24-(L^5*s^5)/120
+(L^6*s^6)/720-(L^7*s^7)/5040+...
結果は
(%o2) [-(4*s^3*L^3-60*s^2*L^2+360*s*L-840)
/(s^4*L^4+16*s^3*L^3+120*s^2*L^2+480*s*L+840)]
と出てきます。結果は同じですね。

> 高次にしないと近似が甘くなることも実感できました。
そうでしょうね。
波形の近似、同じ近似精度なら、できるだけ項数を少ない近似法が重要ですね。

> 5次以上の高次の解を求めるには、何かテクニックがあるのでしょうか。
5次以上の方程式は一般的な解法がないことは証明されているようです。
しかし、解の公式はなくても、
sの実係数多項式のゼロ点は、「実数と共役複素数のセットから構成される」ことが証明されているため、これを利用して現実的な極やゼロ点を分離する理論が確立しているかと思います。
色々なフィルター設計論や回路合成論や自動制御の伝達関数設計論で複雑なsの複素関数合成論が研究尽くされていて、極やゼロ点の最適配置法や部分分数展開法(厳密法法や数値計算による近似法)などの専門書の名著などが出ていますね。もう古典論になっていますので研究テーマとしては取り上げられる事もなくなり、専門書は絶版になっている可能性大ですね。(大学や大学院の単位制の講義には内容が多くて、教科書としてはなじまない為、売れない名著は絶版になる運命にあるのでしょうね。内容のない15回の講義で終わる教科書が多いですね。)
洋書の分厚い古典制御論や回路合成論(複素関数論)などの名著で詳しく扱われていたようですね。今は何でも計算機でリアルタイムの時間領域で制御する現代制御論が主流になってしまった嫌いがありますね。

昔は手計算で計算していたテーラー展開やフーリエ級数展開やパデ関数もパソコンで簡単にあっという間に計算できるようになり、3変数以上の連立方程式はもはや手計算で解くことは無くなりました。

投稿日時 - 2008-05-02 22:15:09

ANo.10

矩形波f(t)のフーリエ級数展開でも基本周波数の10~20倍以上の項まで加え合わせないと矩形波らしくなりません。今回はデルタ関数ですから、さらに急峻な変化を伴う波形ですので、加え合わせる項数を増加しないと急峻なデルタ関数を近似出来ません。パデ近似ではフーリエ級数より少ない項数(次数)の整式の有理関数で近似出来ますが、それでも分子、分母の次数を上げないとデルタ関数のような急峻な波形の再生は難しいと思います。パデ近似はマクローリン展開のn項の和と一致する関数でパデ近似の分子と分母の次数の和がnとなる関係にあります。

また、パデ近似の関数も分子の次数が分母の次数より低い近似関数を使うことで、部分分数展開で定数項が発生しなく出来ますので、そういったパデ関数を使ってみるのも一策かと思います。
いずれにしてもラプラス逆変換で時間領域の矩形波が再現できるわけですから、パデ関数の分子、分母の次数をどんどん上げていけば時間領域の波形も再現できるはずです。実用上どこまで分子分母の次数を下げられるかが問題ですね。

「分子の次数>分母の次数」 のパデ関数(分母の次数が3~6次のもの)を書いておきますので
試してみてください。
[-(6*s*L-24)/(s^3*L^3+6*s^2*L^2+18*s*L+24)]
[(3*s^2*L^2-24*s*L+60)/(s^3*L^3+9*s^2*L^2+36*s*L+60)]
[-(24*s*L-120)/(s^4*L^4+8*s^3*L^3+36*s^2*L^2+96*s*L+120)]
[(12*s^2*L^2-120*s*L+360)/(s^4*L^4+12*s^3*L^3+72*s^2*L^2+240*s*L+360)]
[-(4*s^3*L^3-60*s^2*L^2+360*s*L-840)/(s^4*L^4+16*s^3*L^3+120*s^2*L^2+480*s*L+840)]
[-(120*s*L-720)/(s^5*L^5+10*s^4*L^4+60*s^3*L^3+240*s^2*L^2+600*s*L+720)]
[(60*s^2*L^2-720*s*L+2520)/(s^5*L^5+15*s^4*L^4+120*s^3*L^3+600*s^2*L^2+1800*s*L+2520)]
[-(20*s^3*L^3-360*s^2*L^2+2520*s*L-6720)/(s^5*L^5+20*s^4*L^4+200*s^3*L^3+1200*s^2*L^2+4200*s*L+6720)]
[(5*s^4*L^4-120*s^3*L^3+1260*s^2*L^2-6720*s*L+15120)/(s^5*L^5+25*s^4*L^4+300*s^3*L^3+2100*s^2*L^2+8400*s*L+15120)]
[-(720*s*L-5040)/(s^6*L^6+12*s^5*L^5+90*s^4*L^4+480*s^3*L^3+1800*s^2*L^2+4320*s*L+5040)]
[(360*s^2*L^2-5040*s*L+20160)/(s^6*L^6+18*s^5*L^5+180*s^4*L^4+1200*s^3*L^3+5400*s^2*L^2+15120*s*L+20160)]
[-(120*s^3*L^3-2520*s^2*L^2+20160*s*L-60480)/(s^6*L^6+24*s^5*L^5+300*s^4*L^4+2400*s^3*L^3+12600*s^2*L^2+40320*s*L+60480)]
[(30*s^4*L^4-840*s^3*L^3+10080*s^2*L^2-60480*s*L+151200)/(s^6*L^6+30*s^5*L^5+450*s^4*L^4+4200*s^3*L^3+25200*s^2*L^2+90720*s*L+151200)]
[-(6*s^5*L^5-210*s^4*L^4+3360*s^3*L^3-30240*s^2*L^2+151200*s*L-332640)/(s^6*L^6+36*s^5*L^5+630*s^4*L^4+6720*s^3*L^3+45360*s^2*L^2+181440*s*L+332640)]

ある程度の幅を持った孤立矩形波の場合の伝達関数のラプラス逆変換のパデ近似や同じインパルス応答のパデ近似の分子、分母の次数を増加して行った時の応答の変化を調べることで、どの位のパルス幅の孤立矩形波(三角波でもよい)やパデ近似の分母の次数をどの位取らないとインパルスが再生できないかを調べてみるといいでしょう。

投稿日時 - 2008-05-01 14:54:41

お礼

ご回答頂きありがとうございます。
しかも6次までのパデ近似、自分はせこせこ連立方程式を計算していましたが、
どうやったらこれだけ多くの計算を短時間に出来たのでしょうか?
何か特別なやり方があるのでしょうか。エクセルでプログラムを既に
組んでらっしゃるのでしょうか。よろしかったらご教示いただきたく
お願いいたします。

ご教示頂いた中で、簡単そうな
[(3*s^2*L^2-24*s*L+60)/(s^3*L^3+9*s^2*L^2+36*s*L+60)]
について解き、L=1,2,5のときについてエクセルで描かせてみました。
振動は大きいですが、L=1,2,5の所に山ができ、一応インパルス応答を
再現できました。ありがとうございます。
高次にしないと近似が甘くなることも実感できました。
ただし、高次、特に5次以上は根を求めるすべが代数学的にないと認識しております。
4次方程式はフェラリの解法というものを使えば出来るそうですが、
小生3次方程式までしか勉強しておりません。
5次以上の高次の解を求めるには、何かテクニックがあるのでしょうか。
ご存知でしたらご教示いただきたくお願いいたします。

投稿日時 - 2008-05-02 19:07:01

ANo.9

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

[1] パデ近似自体が全然駄目ということではないでしょう。
 インパルス応答のうちでδ(t)を無視した波形を見れば、δ(t-T)の(なだらかだし、振動しているけど)近似になっている筈です。  
 ANo.6のお礼で、

> δ関数を次式で近似しました。

と仰るのは、G(s)の逆ラプラス変換をガウス関数で平滑化したものを計算してみたということでしょう。(本来は平滑化をしようがしまいが関係ありませんが、δ(t)じゃexcelに乗らないからこうなさったのでしょうね。)
 その結果、一番高いピークはt=0に出るにしても、t=Lの所にもそれなりに盛り上がりが観察された筈です。ガウス関数の幅を広げて平滑化を強くしてやるにつれて、t=0のピークがどんどん下がって目立たなくなる一方で、t=Lの盛り上がりはあまり低くならない筈です。
 つまり、「t=Lの盛り上がりの幅が比較的小さくて、しかも(山の高さではなく)カーブ下面積で見てt=0の山に比べて大きくなっていてくれさえすれば、急な変化をしない入力に対してならば、G(s)は無駄時間要素の代用にならなくもない」というのが、パデ近似の信号空間に於ける意味でしょうね。

[2] | s|→∞でG(s)→Kのとき、ランプ応答ならば(s^(-2)を掛け算するわけだから)、Kδ(t)の成分は信号空間ではKtとして現れる。Kが他の成分に比べて相対的に小さければ余り目立たないでしょう。でもインパルス応答だとKδ(t)がモロに目立つので、高周波成分を含む激しく変動する信号には使えない。実用を考えると、せめてステップ応答で使いものにならないと面白くないんじゃなかなあ。
 性能を上げるために、近似になってるsの範囲を広げれば、Kδ(t)に比べてt=Tの山が大きくなり、しかも細くなって、δ(t-T)に近づく。しかし、そのためにはパデ近似の次数を高くする必要があり、要素Gを構成するのが大変になる。そういう構造の問題です。

[3] もしご質問が実務の話だったら、パデ近似にこだわらないという選択もあるはず。
 たとえばANo.8のやり方ですと、振動は残るがδ(t)は完全に消えます。
 ANo.5の近似は、δ(t)が現れないし、次数nが案外小さい段階でかなりδ(t-T)に似ます。それに、全く振動しません。山の幅がそのまま、カットオフされる高周波の閾値を表している訳ですから、nを幾らにするかの設計が簡単。総合的に、パデ近似より良いと思います。弱点は無駄時間が入力信号の周波数によって少々違うことですけど、パデ近似だって波形が崩れるのは同じこと。
 いや、そもそも無駄時間要素は、電気信号なら例えば表面弾性波を使って素直に実現する方法もあるから、何が何でも近似をするしかないという訳ではないんだし。

投稿日時 - 2008-05-01 14:31:21

お礼

ご回答頂きありがとうございます。
現在は実務ではなく、将来の向学のために制御工学を勉強中です。
近い将来必要になることを見越してです。
どうせ勉強するなら実務で使えるものでないと意味がないと考え、
実務想定で勉強中です。
>ANo.8のやり方ですと、振動は残るがδ(t)は完全に消えます。
参考になります。ANo.6の偶数次の近似式と、ANo.2お礼にある
奇数次の近似式とを足して2で割るということですよね。
試してみます。
>ANo.5の近似は、δ(t)が現れないし、次数nが案外小さい段階でかなりδ(t-T)に似ます。...
せっかくご教示いただいたのですが、高等すぎてよく分かりません。
ご面倒でなければ、ANo.5で省かれていた証明をご教示いただければ幸いです。

投稿日時 - 2008-05-02 18:55:25

ANo.8

ANo.7を書き込んだ直後に思い至ったのだけれども、ANo.6の偶数次の近似式と、ANo.2お礼にある奇数次の近似式とを足して2で割れば、|s|→∞のときに0になる近似式になる。当たり前だあ。これでも「パデ近似を使ってる」ことには違いない。展開したらどうなるでしょうかね。

投稿日時 - 2008-04-30 17:39:14

ANo.7

ANo.4の追記です。

 |s|→∞のとき、もしG(s)が0以外の定数Kに収束するんだと、逆ラプラス変換で(δ(t-T)ではなくて)Kδ(t)が現れる。ってことはつまり、G(s)が有理関数の場合、分子の次数は分母の次数より低くなくてはならない。
 なので、パデ近似のアプローチでは、逆ラプラス変換したときにδ(t-T)を近似するような式を作るのは無理だと思う。

投稿日時 - 2008-04-30 17:20:41

お礼

ご教示ありがとうございます。
ANo.6お礼欄に記載させていただきましたが、やはり
δ(t-T)ではなくてδ(t)が現れるため、求めようとした結果になりませんでした。
ステップ応答やランプ応答は求めようとした結果が出るのに
なぜインパルス応答だけは、得ようとした結果にならないのでしょうか。
これでは近似の意味がないですよね。
しかし、制御工学の教科書には無駄時間要素は有理関数のパデ近似する
と書いてあります。そもそもインパルスだけは必要がないのでしょうか。

投稿日時 - 2008-05-01 11:47:25

ANo.6

#2です。
>無駄時間要素をe^(-Ls)とすると3次/3次のパデ近似
> e^(-Ls)
> ={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}
> /{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120}
この分子,分母のsの最高次係数の比が「-1」になり、部分分数展開の定数項が「-1」となって-δ(t-L)が発生します。
これはs→∞でバデ近似関数→-1になることに起因しています。つまり、
これがt=Lで「-δ(t-L)」の発生要因になります。
一方、s→→∞でe^(-Ls)→0+ になります。つまり、
これがt=Lで「δ(t-L)」の発生要因になります。

この矛盾をなくするには、パデ近似関数の分子・分母の次数が偶数次の展開を使用する必要があります。つまり、e^(-Ls)のパデ近似関数として
(2次)/(2次)の
(s^2*L^2-6*s*L+12)/(s^2*L^2+6*s*L+12)
または
(4次)/(4次)の
(s^4*L^4-20*s^3*L^3+180*s^2*L^2-840*s*L+1680)
/(s^4*L^4+20*s^3*L^3+180*s^2*L^2+840*s*L+1680)
関数を用いれば、
分子,分母のsの最高次係数の比が「1」になり、部分分数展開の定数項が「1」となってδ(t-L)が発生させることができます。
これらのパデ関数ではs→∞でパデ関数→1>0になり、このことがt=Lでδ(t-L)を発生することを意味します。

ということで無駄時間要素(遅延要素)のパデ関数として
(sの偶数次)/(sの偶数次)の近似関数を使う必要があると言うことですね。
そうすると、近似関数を部分分数展開は定数項として「1」を分離すれば、残りは普通の部分分数展開ができるようになり、ラプラス逆変換も可能になるかと思います。

投稿日時 - 2008-04-30 02:42:48

お礼

ご回答頂きありがとうございます。
ご教示頂いたe^(-Ls)のパデ近似関数として
(2次)/(2次)の
(s^2*L^2-6*s*L+12)/(s^2*L^2+6*s*L+12)
でインパルス応答を求めてみました。
解は共役複素数の2点が求まり、やはり複素平面上では、
分母(極)と分子(零点)は左右対称でした。
一般解をσ±jωとすると
G(s)={(s+σ+jω)(s+σ-jω)}/{(s-σ-jω)(s-σ+jω)}
 ={(s+σ)^2+ω^2}/{(s-σ)^2+ω^2}
とおけるので
G(s)=A+{B(s+σ)+Cω}/{(s-σ)^2+ω^2}
とおいて部分分数展開し連立方程式を求めると
A=1
B=4σ
C=4σ^2/ω
となりました。
これをL=1,2,5の場合について解きδ関数を次式で近似しました。
δ(t)=lim{n→∞}[(n/π)^0.5]EXP(-nt^2)
エクセルで描いてみたところ、δ(t-L)なので
ピークはL=1,2,5で出るはずが、
すべてt=0のところになってしまいました。
無駄時間要素においてパデ近似したもののインパルス応答は求められないのでしょうか。

投稿日時 - 2008-05-01 11:39:51

ANo.5

 無駄時間要素を有理関数で近似するという発想が面白いんで、ちょこっといじくってみました。

γ(n,T,t) = (t^(n-1)) exp((1-n)t/T)
という、t≧0で定義された実関数を考えます。ここに、T≧0、nは正の自然数です。
 すると、γは連続関数であって、
・ γ(n,T,0)=0
・ t>0のときγ(n,T,0)>0
・ t→∞のとき、γ(n,T,0)→0
・ t=Tに唯一の極大を持つ。
・ t=T±T/√(n-1)に変曲点を持つ。
という性質があります。
 また、γのグラフの面積
A(n,T) = ∫γ(n,T,t)dt (積分範囲は0~∞)
は当然ながらnとTだけの式で書けます。(Γ関数を使えば具体的に計算できますが、めんどくさい(ってか、実は間違えそうな)のでパス。)

 そこで、γをAで規格化したもの
f(n,T,t) = γ(n,T,t)/A(n,T)
を考えると、n→∞のとき、f(n,T,t)→δ(t-T) って寸法でさあ。(え?証明? えーと、めんどくさいのでパス。<またかよ)この関数、結構性格が良さそうだってことは、幾つかのnについてグラフを描いてみればお分かり戴けるのではないかと思うんだけどな。

 一方、
G(n,a,s) = 1/(s+a)^n
の逆ラプラス変換は
g(n,a,t) = (t^(n-1)) exp(-at) / (n-1)!
だから、
γ(n,T,t) = (n-1)! g(n,(n-1)/T,t)
ゆえに
f(n,T,t) = ((n-1)! / A(n,T)) g(n,(n-1)/T,t)
つまり、fのラプラス変換は
F(n,T,s) = ((n-1)! /A(n,T)) / (s+(n-1)/T)^n
である。

 という訳で、無駄時間要素 exp(-Ts)(ANo.4では「遅延」と書きましたが、「無駄時間要素」と呼ぶ方が良いですね。)は F(n,T,s) (n→∞)で近似できる。つまり、単に1次遅れ要素 1/(s+a) をいっぱい直列してやれば、(従って、exp(-at)をいっぱい畳み込めば)T=(n-1)/aの無駄時間要素に似てくる。こりゃ知りませんでした。

投稿日時 - 2008-04-29 23:39:56

ANo.4

 ANo.2に付けられたコメントについて、またしても横入り御免なすってのstomachmanです。

 お示しの近似式は|s|→∞のとき-1に収束するでしょ。つまりG(s)は|s|がごく小さい時を除いて
G(s)=-1
である。これが-δ(t)が現れる理由です

 G(s)で遅延要素を近似しようとするなら、「近似として有効な範囲を越えたら0に近づく」ような仕掛けを入れる必要があるようです(が、そこまでやったら有理関数ではなくなって、m近似する意味が失われる?)。

投稿日時 - 2008-04-29 16:42:58

ANo.3

 横入り御免なすって。えとですね、デルタ関数が出て来てもちっともおかしくありません。物理的にも、ごく普通にあることです。というのは、「入力を、なんにもしないでそのまんま出力する系」の伝達関数 (modulation transfer function)こそがδ(t)だからです。
 なので、-δ(t)とは「入力の符号を逆にして出力する系」の伝達関数です。
 ANo.2にある部分分数分解(チェックしてませんが)は「3つの項がそれぞれ伝達関数を表していて、これら3つの伝達系が並列に繋がって出来ているのがこの系である」と解釈される。何の不思議もありません。

 一般に、ラプラス変換で表した伝達関数G(s)の分子の次数は、その分母の次数より低い場合もあれば、同じの場合も、高い場合もあります。伝達関数が積分的な作用をする(例えばコンデンサ)なら、分子の方が分母より次数が低い。微分的な作用をする(例えば微分フィルタ)なら、分子の方が分母より次数が高く、この場合は逆ラプラス変換をするとδ関数の微分(微分演算子。これを関数f(t)と畳み込み積分すると、df/dtになる)が出てきたりします。
 そして、一定時間の遅延を掛けるとか、位相を変えるだけの作用をする伝達関数は、分子と分母の次数が同じになる。例えばe^(-Ts)なんてのは次数0で、逆ラプラス変換すればδ(t-T)、つまり遅延を表す伝達関数です。

投稿日時 - 2008-04-29 04:36:52

ANo.2

数式的には
> G(s)=-(s-a){(s-σ)^2+ω^2}/[(s+a){(s+σ)^2+ω^2}]
> =A/(s+a)+{C(s+σ)+Dω}/{(s+σ)^2+ω^2}
=-1+{A/(s+a)}+{C(s+σ)+Dω}/{(s+σ)^2+ω^2}
と部分分数展開されます。

ラプラス変換
∫[0→∞] x(t)e^(-st)dt=-1
となるx(t)は
x(t)=-δ(t)

このδ関数(デルタ関数)は
[Diracのδ関数、インパルス関数]
http://ja.wikipedia.org/wiki/%E3%83%87%E3%82%A3%E3%83%A9%E3%83%83%E3%82%AF%E3%81%AE%E3%83%87%E3%83%AB%E3%82%BF%E9%96%A2%E6%95%B0
と呼ばれる関数です。
物理現象として無限大の振幅の電圧や電流などの物理量は存在しませんから、それをラプラス変換したs領域で定数項が現れることはありません。
つまりラプラス逆変換の対象となるG(s)の分子は分母よりsの次数が低くてはいけないということです。
あえてG(s)の分子と分母が同じ次数としたため、部分分数展開で定数項が発生し、時間t領域にデルタ関数(∞振幅、幅ゼロの関数、面積1の関数)がラプラス逆変換で発生するのです。

投稿日時 - 2008-04-28 18:35:16

お礼

ご回答を頂きありがとうございます。
部分分数展開で定数項を入れることをご教示頂きありがとうございます。
おさっしのことかもしれませんが、小生の質問は、無駄時間要素
をパデ近似(3次/3次)したもののインパルス応答を求めようとしています。

無駄時間要素をe^(-Ls)とすると3次/3次のパデ近似は次のようになります。
e^(-Ls)={1-(Ls)/2+(L^2*s^2)/10-(L^3*s^3)/120}/{1+(Ls)/2+(L^2*s^2)/10+(L^3*s^3)/120}
上式の分母分子の3次式をそれぞれカルダノ法でとくと1つの実根と
2つの虚根(共役複素数)が得られます。
これらを、複素平面上で表わすと分母(極)は左半面、
分子(零点)は右半面で分母分子は左右対称であることから
一般解を実根a、複素根をσ±jωとおき、分子多項式の係数の正負に
注意すると伝達関数は次のようにおけます。
G(s)=-(s-a)(s-σ-jω)(s-σ+jω)/(s+a)(s+σ-jω)(s+σ+jω)
インパルス応答g(t)はg(t)=L^-1[G(s)]なので上式の
ラプラス逆変換を求めようとしています。

無駄時間要素のインパルス応答は近似を使わないで求めると
L[f(t-L)]=F(s)e^(-Ls)を利用してg(t)=δ(t-L)となります。
これは、無駄時間Lだけシフトしたδ関数になります。
しかしながら、部分分数展開の過程で-1の項がでるということは、
おっしゃるとおり-δ(t)の項が出てきてしまい。正負反転してしまいます。
小生の求め方にどこか間違った点があるのでしょうか。
ご指摘、ご指導いただければ幸いです。

投稿日時 - 2008-04-29 14:13:25

ANo.1

分子の次数が分母の次数より小さくなるように, 先に分子を分母で割っておく.
でも, 1 を逆変換すると δ(t) になるんじゃなかったっけ? これで意味のある答えになるのかなぁ.

投稿日時 - 2008-04-28 18:29:09