重みつき残差法による分布定数系のモデリングに関 する研究

重みつき残差法による分布定数系のモデリングに関 する研究

今井, 純

九州大学工学研究科電気工学専攻

https://doi.org/10.11501/3088157

出版情報：Kyushu University, 1991, 博士（工学）, 課程博士 バージョン：

権利関係：

第4章 3次スプライン関数による近似モデル

4. ₁ 緒言

重みつき残差法では近似の目的に応じた座標関数の選択が重要である。

区分的多項式の一種であるスプライン関数は関数近似において一般的に用 いられており 特にl次スプライン関数を用いる有限要素法は偏微分方程 式の数値的な近似解法として普及している。 しかしながらこのような低次 のスプライン関数を用いると近似モデルそのものの次数は高く なって制御 系設計や推定問題などへの応用には不都合な場合がある。 そのため， 例え

ば直交多項式， 固有関数や高次のスプライン関数など， より滑らかな座標 関数がそれらの目的には用いられてきている。 特に3次スプライン関数の 様々な推定問題への適用に関する研究は活発に行われている14)ー18)37)。

重みつき残差法の近似モデルの行列の要素の計算は， 数値積分により処 理される例が多い。しかしパラメータ推定や適応制御などのようにパラメー タを変更しながらモデルの構成を逐次繰り返さねばならないような場合に は必要な計算量が膨大になるという問題がある。 こうした問題を回避する ためには問題の行列の計算を数値積分に任せるのではなく， 事前に解析的 に求めておく ことが特に有効である。 しかしながら高次のスプライン関数 はその表現が複雑になるためそのようなモデルの構成を具体的に示した例 は見かけないようである。

最近， 有限要素法における行列の計算に数式処理を用いて記号的に積分 を行わせ， その後数値計算させることで計算の精度の改善や負荷の軽減を

図る試みが行われている。

本章では1次元分布定数系の3次スプライン近似において数式処理を援 用した近似モデルの行列の解析的な構成を提案する。 これによりパラメー

タ推定問題のように数値計算のみによる処理では膨大な計算量が要求され る分野にもスプライン関数近似が適用できるようになる。デイリクレ型境界 入力を持つ放物型l次元分布定数系の3次スプライン関数による近似モデ ルが一般によく用いられている固有関数近似と比較して優れた精度を有す ることを時間領域ならびに周波数領域における誤差評価の比較により示す。

4. 2 対象システム

本章で対象とするl次元の空間領域をr2== (0， e)とする。 。上で定義さ れた状態、u(t，x)，xεr2，t>0 の時間的空間的挙動が、次式で与えられる放物 型分布定数システムで記述されるとする:

Z

M=AUM十川仇 t>山r2， (4.1) fu(t，x) == f(t) ， t> 0

u(O ， x) == uo(x)， xεQ ここでAは楕円型偏微分作用素、 rは境界作用素で

A _{ψ :戸=}α仰刷(μωZ刈)

2宗P +吋γ村仰(μ暗Z

I 川 |い

F伽州凶川川( 仁川0的)ト一べ (1ト一

F仇ωLρ)祭(t附川，0的) I

fu川叫仰(tい川，xめ)⁼

^l |

^f_rR ^R

^J |

u ^� := ^'-

I !

ßRU(t，e)+ ( トグR)

髭

^{( υ )}

| I

( 4.2) ( 4.3)

でありo ::; ßL ::; 1 ， 0::; ßR ::; 1 ，α(x) > 0 ， c(x)三0である。境界条件(4.2)式 右辺に現れている境界入力了(t)== [fL(t) fR(t)]Tは超関数的な入力影響関数

を持つ分布入力へ変換することができるのでOと仮定する。

スプライン関数は区分的な多項式として表現され、座標関数として用い るには基底の修正が必要なので近似モデルの計算は固有関数の場合に比較 しでかなり面倒になる。以下では同次境界条件を満たす座標関数の構成と 近似モデルの行列要素の解析的構成について述べる。

4. 3 座標関数の構成

k次 スプライン関数とはQ上に配置された節点で区切られる各小区間で k次多項式でかつ節点上でk - 1次までの導関数が連続であるようなもの をいう 11)。以下では領域n ⁼⁼ [0， f]上に 境界も含めてN個の節点をo ⁼⁼ Xo <

Xl < ・・・ くXNー^{1 ==} f， Xi ⁼⁼ ihとする。ここで hは節点の間隔f/(N -1)であ

り、これらの節点はQをN - 1個の等しい長さの小区間へ分割する。この ような節点上で定義されたスプライン関数全体か らなる部分空間の基底と しては通常Bスプライン関数が用いられる。等間隔の節点の場合、 正規化

された3次のBスプライン関数は次式のように 表される 38).

BJ(x) ^:==

(X - Xj_2)3 6fh3

h3 + 3h2(x-Xj_l) + J 3h(x-xjー1?J - 3(X-Xj_l)3 J ， zε[Xjー1，Xj) 6fh3

h3 + 3h2(Xj+l-X) + 3h(Xj+l-X)2 - 3(Xj+l-X)3

6fh3 . - ， Xε[Xj， Xj+I) zε[Xj-2' Xj-I)

zε[Xj+l， Xj+2)

Z ダ[Xj-2'Xj+2) ( 4.4 ) Bj(x)はFig. 4.1に示すように 、節点zjで頂点となるような釣り鐘型のグラ

フとなる。 4つの小区間[Xj-2'Xj-l]，. . . ， [Xj+l， Xj+2]上でそれぞれ3次多項式 で表され、それ以外では恒等的に0である。また、Bj(x)は必ず0もしくは

正の値を取り、その積分j

二

Bj(x)dxは節点間隔 hに 等しい。

(Xj+2 - X)3 6fh3

ハU

考えている空間領域n ⁼⁼ [0刈の上でOでない値を取るBスプライン関数は Bj(x)，j ⁼⁼ -1，0，..けNのN+2個である凶)⁰ Bj(x)のうち、j⁼⁼ 2， 3，・・・ ，N -

2については境界XO，XNー1上でその微係数も含めて0となるので問題はない が、Bj(xo)，j⁼⁼ -1，0，1， Bj(XN-I)，j ⁼⁼ N - 2， N - 1， Nはそうで ないため同 次境界条件を必ず しも満足せず 、これらをそのまま座標関数として用いる ことはできない。そこでB-l(X)，BO.(x)， Bl(X)およびBN(x)， BN-l(X)， BN-2(X) 注2)境界条件を適切に処理するため、ここに限ってQの外の節点Xj^:== jh， j ⁼⁼

-2， -1，N，N +1を考えている。

兄!て 9/1 0 (X)[H

トミ+

. .，..圃・3

K

F吋+

J3

J2

�

J3

トミ

応、

同

ま議区

，，，\

ヘャ II'\

Ì'\

tく

手《凶

ゲコ

ヰと。

外j

話也;

活きヰさ

F→

寸4 b.O

�

をそれぞれ 線形結合 することによ・り 同次境界条件を満たす2個ずつの座標 関数仇(x)，仇(x)および<ÞN-l(X)，ゆN(X)を作り、次に示す ように考える39)。

αLBO(x) + bLB-l(X)， i == 1 cLB-1(x) + dLBl(X)， i == 2

仇(X):== < B←l(X)， 3三t三N-2 (4.5)

cRBN(x) + dRBN-2(X)， i == N-1 αRBN-l(X) + bRBN(x)， i == N

こ こで係数α*，b*， c*，ム(ただし、* == L， _R ) は向次境界条件r仇==0と、正 規化の条件(ん仇dx == h ==一定)より定まる線形方程式を解くことにより次

のように求められる。

α輩 :== 6 ß*h+3(1-ß*) 皐 2ß*h+ 9(1 -ß*) 払:==-2 4 ムh

事 2ß*h+ 9(1 -ß*) ムh -3(1 -ß*) ら:==-12

11ムh+ 36(1 - ß*) ß*h + 3(1ーム) ι:== 12

原 11ムh+ 36(1 -ß*)

( 4.6)

ただし、 N三4である。このようにして定まっ たゆi(X)は常に同次境界条件

を満足するようになるが、その ため境界条件のパラメータ(ßL，ßR)に依存 し たものになっている。 仇(x)の様子をグL == ßR == 1， N == 8の場合について

Fig. 4.2に示す。

キ弘、

υつ

8 0

Xo Xl X2 X3 X4 Xs X6 X7

X

Fig.4.2 ßL ⁼⁼ ßR ⁼⁼ 1について3次Bスプライン関数を修正した座標関数( ^{N ==} 8 ) ^σコ

心2

4. 4 行列要素の解析的構成

パラメータ推定問題のようにモデルの構成を繰り返すような場合には (2.56)式に現れる行列M，

伊

^，A<ÞT

)

，(φ?ψ) はパラメータを含んだままの解析 的表現で求めることが望ましい。しかし、 3次スプライン関数は表現が複 雑になるため、手計算だけでそれを求めるのは困難になる。以下では数 式 処理の援用を前提とした、そのような行列要素の計算の方法について説明 する注3)。

まず(ゆi，Acþj):==

ρ

i(x)Aめ(x)dx という計算を実際にはどのような tと jの組み合わせについて求める必要があるかについて考える。 Aが自己共 役すなわち(札Aゆj)== ( A仇?ゆ'j) ならばi，jに関する対称性が成り立つこと を利用して、 t三 jなる成分のみを計算する。またli - jl三4 のとき被積分 関数 仇(x)AゆIj(x)は恒等的に0になるのでli - jl _� 3 の場合のみ (仇，AゆIj)を 求めればよく、Fig. 4.3に示すように行列

何

^，AφT

)

は7重対角行列になる。

スプライン関数を用いている場合の要素計算は次のように考えれば良い。

すなわち各小区間についてそれぞれ被積分関数の陽表現を求め、不定積分 を計算したのち積分区間の代入により定積分を得て 各小区間についての総 和を取れば結果が得られることになる。行列要素 (仇，Aゆ'j)の計算において は、( 4.5)式のような境界条件による基底の修正が 行われているため、境界 が関係しない3三t壬N-2かつ3三j_� N-2の場合とそれ以外の境界が関 係する場合とに分けて考えなければならないo A:== 1 (恒等作用素)の場合、

伶

^JφT

)

はグラミアンM になるが、ここでは実例としてこの場合を説明す る。 Z == Jのときは被積分関数が零でない値をとる 4つの小区間[Xi-2' Xi-l] ，

注3) ここに述べる行列要素の計算だけではなく、( 4.4)式のBスプライン関数 40)や(4. 6)式の計算などにも数 式処理を効果的に適用できる。

∞ロoz一℃口oub弓503 ロ'Hω∞ロ

oυ 最長室一語e73〉誌に

め .

叩.∞一』

ミミ

寸 一寸 一」

一寸

_.

^-i

_._.

^�

一一

^L

_.

^JI

4

ミ

ーし| 円lH〈NJ之のJ〈マはく

.ぐh、

11

ば〉

寸4

止

寸4

止 止

寸4 ばコ N σつ

•

<:'> F→

.小、

[Xi-l' Xi]， [Xi， Xi+l] ， [Xi+l' Xi+2]についてのみ各々積分を考えれば良い。

川:=

[

_{(ダーl(X))2 dx}

=

[��' (

^{�X - Xi-3)}

^� r

_ \ ^6eh3 +

[ � �'C3て了'C3C3が山山3斗愉+刊州3幼h件2句(

一 ωt 一 )

^{dx (4.7)}

十日

h3 + 3h2(XiーX) + 3h仇- x)2- 3(Xi- X)_6eh3

� r

J +

[H' (川 -

^x)

^� r

\ ^6eh3 J

またj==i-l のとき考えるべき小区間は3区間でよい。

(<Þi，<Þi-l):=

II

Bi-1(X)Bi-2(X)dx rXi-2 (X - Xi_3)3

ムー3 ^6eh3

h3 +3h2(x-Xi_3)+3h(x-Xi_3)2ー3(x-x←3)3dz 6eh3

1:1山印一勾九山一-2

6eh3

h3 + 3h2(Xi_l - X) + 3h(Xi-l - x)2 - 3(勾ー1- x)3dz 6eh3

+

(

^X

‘

h3 + 3h2(Xi - X)ギ制Xi- X)2ーめi- x)3

L・ 1 ^6eh3

(Z川- ^X)3Jー 6eh3 ^��

(4.8)

j == i -2， i - 3の場合も同様である。このような多項式の積分の操作は数 式処理に向いており、計算機上で容易に処理が実現できる。Xi :== ih， Xi-l :==

(i -l)hなどを代入すると以上の計算結果は(4.9.a)式になる。

Mi，j ==

151h/315，

39 7 h/168 0，

h/42，

h/504 0，

O

z-j==Oのとき li -jl == 1のとき li -jl == 2のとき li -jl == 3のとき li -jl三4のとき

(4.9.α)

ただし、Mi，j は行列Mのt行j列の要素である。

tまたはjが 1， 2，N-1，Nの場合 は、境界が関係してくるため境界パラ

メータグL，ßRを含むことになるが、 同様の手順で各要素について行うと次 のような結果になる。

ー ー 1 2114 () 43 ^{1 _} () -1 M1.1 ^:== I一一7:a1.+一一似仇+^�=�b

�

_{I h}

1 8820 --ú . 840 --^{... r} U • 252 んl

1 43 354 1 ^_ 1 ^{_ _} 1

Ml，2=M1 1:=|一一一αLCL+一一一αLdL+ ^":，， bLcL + ^�-， bLdL

I

_h

1 1680--u-u . 1680--u--u . 252-u-u . 84-u 心 | M^1._v ^{3 == M3}_v，" .1 ^:==

I 土

_aL+ _�1，� bL

I

_h

1 42 - U • 5040 -u 1

M1，4 ⁼⁼ M4，1 ^:==

I

15040 _r::

^}

_{-! {\aL} -u

I

1 _h (4.9.b)

_ _ 1 1 ^{'l 1} 4193 ^{_'l 1}

M2.2 ^:== I一一c1.1252 -ú・42-+ ^，-� cLdL U--U・8820--ú +一一-dlh¹

1 1 397 ^_ 1

M2.3 ^� == M3.2 v，� ^:== I一一一CL+一一1 5040 � ^. 1680 ..� dLI h ^{� 1} M2.4 ⁼⁼ M4.2 _ι，� ^:==

I

�1. dL

I

_h

124 � 1 M2.5 _υ^{== M}_v，" 5^.2: ⁼⁼

1

_{1 5040 _�1.� dL}

I

_h

--U 1

一般のAの場合、Fig. 4.3に示すようにN三7 の場合、 3�i� N-2， 3�j三N-2のときは、境界が関係しないのでNおよび丸jを含んだ一般 式を導くことができる注4)。 そのようにして ， Aが自己共役のとき20通り、

自己共役でないとき35通り の要素の計算で十分である注5)。以上のような結 果を導いておくことで、 システムのAゃr， bjなどに表れるパラメータを 単純に 代入することにより近似モデルの構成に必要な行列が直ちに得られ ることになる。こうした代入計算や(2.56)式の逆行列計算など は数値的に容

易に実行される。

注4) N三6について は全てのli-jl三3なる要素を求める必要がある。

注5)節点が等間隔でない場合 は全要素についての計算が必要となり、 その個 数は自己共役のとき 4N-6 ， 自己共役でないとき 7N-12 になる。

4. 5 数値例

本節ではデイリクレ境界入力をもっ系について誤差の空間分布の様子と その最大絶対値の周波数特性を示し、試行関数としての固有関数と3次ス プライン関数による近似の性能の違い を明らかにする。

( 4.1)一(4.3 )式 において、領域は1次元でn = (0，1)とし、 A=θ2/θd，ru = [u(t，O) u(t， 1)]' (デイリクレ型境界条件)と する。

θu θ2U

-- - 一τ+b( x ) f ( t )， 0 < x < 1， t > 0 θt 8x2

u(t，O) = u(t， 1) = 0， u(O， x) = 0

(4.10)

いまó(x) をデイラックのデルタ関数としてb(x)=グ(x-1)の場合を考える。

このような入力は u(t，x)の空間分布 における x= 1-0 において-f(t)の不 連続を生じさせる働きをすると 考えると、 u(t，l)=O より、limx-+lーQu(t，x)=

f(t)である。するとこの問題は次のような境界値問題と等価である。

θu δ2U

一一=一τ ， 0 < xく1，t_> 0 δt 8x2

u(t，O) = O，u(t， 1) = f(t)， u(O，x) = 0

がり)+グ(x - 1)=0，ゆ(0)=ゆ(1) .= 0の解を次のよう に考える:

(x， 0くx<l 0(z )=

{

^{;? ;二1}

(4.11)

(4.12)

座標関数として(4.6)式の境界条件を修正したBスプライン関数を用い て 近似モデルを構成する。 M，

(

^φ，AφT

)

，(φ?ゆ)を3 2ピットパーソナルコン ピュー タ上で動作する汎用数式処理システムREDUCE 3.3で求め、FORTRAN コードの形で 生成させた結果を用い 数値計算を行う。この数式処理を行う ためのソースプログラムのリストを付録. Bに示 す。

3次スプライン関数を座標関数として 用いた場合と 固有関数{sin n7rX， n = 1， .・・ ，N}を用いた場合と について、 それぞれん(t) := sin 97rt ，N = 6，7，8と

してシミュレーションを行った。 シミュレーション結果をそれぞれ厳密解と 比較し、空間についての最大誤差maxxED lu(t，x) - uN(t，x)1を各時刻で求め たものをFig. 4.4に示す。

加法的誤差の空間的最大絶対値ムfriax(ω) :== maxxED IG(jω，x) - GN(jω，x)1 の周波数領域におけるプロット(近似次数N == 4 ，8，16 )を3次スプライン近 似と固有関数近似について それぞれ Fig. 4.5の実線と破線に示す。低周波域 へ向かうにつれ誤差が減少していく割合は固有関数近似の場合、20dB/dec であるが、3次スプライン関数近似の場合、40dB/decと大きい。 また十分 低いωについては固有関数近似の場合、近似次数を2倍にしたとき誤差は

1/4になるが、3次スプライン関数近似の場合、1/20になっている。

以上のデイリクレ型境界入力の例においてはスプライン近似の性能が優 れていたが近似の性能は入力影響関数の形に依存し、必ずしも著しい性能 の向上が見られるとは限らないが、点入力や一様入力など著者らが調べた 限りでは固有関数の場合と少なくとも同等かそれ以上の精度が得られて いる。

o

T-'I

o

://

/

℃HZ

/ /

ノ/

何削州柑SE世(も州制山総坦友裂け帆雌SE山M叩.叩・∞一向

+ー

\

\

" "

" "-

" "-

\、 \、

\、 、\

、、句、\ ー\、

、、、、 \、\ 、

、、、-- 、 、、、 \、"

『、、』 、、 、\

、、 、h旬 、、、、 、\、

----二一ミ\_{�-、子、〉、、}

九 ミ ミミ

\

。HHZ

。ロ

.""，吋 F田吋巳4

� α3

ro ü

咋「 ・�

6ぞ Sぢ

o c c c

o

3XXBlli

tn o o o

I (X'l)NU-(X'l)n

�

NIC)[

�

寸

o

I (X'

^ro

D ^v I

^03xXBlU

v)

o �

てず0¹

�尽

�

1回111区〈 ミ2土日

3

_4く^{千E二コミ}

同長ばつ 寸4 '1"'"'-1

K b 4

。 O

�

H IC-

4. 6 結論

l次元分布定数系の3次スプライン近似において数式処理を援用した近 似モデルの行列の解析的な構成を提案した。 これによりパラメータ推定問 題のように数値計算のみによる処理では膨大な計算量が要求される分野に もスプライン関数近似が適用できるようになる。 デイリクレ型境界入力を 持つ放物型l次元分布定数系の3次スプライン関数による近似モデルが一 般によく用いられている固有関数近似と比較して優れた精度を有すること

を時間領域ならびに周波数領域における誤差評価の比較により示した。

第5章 分布定数系のパラメータ同定問題への応用

5. 1 緒言

あるモデルを用いて制御系を設計するためには、仮定した数学モデルが 現実のシステムに最もよく適合するようにパラメータを決定する、 いわゆ るパラメータ同定の手続きが必要である13)41)-43)。分布定数システムに対 する同定問題についてはこれまで多く研究され論じられてきたが4 4)45)、そ のほとんどは偏微分方程式に現れるパラメータの決定を取 り扱ったもので ある。有限領域で定義された分布定数システムの挙動においては、領域内 部のみならずその境界における特性も重要であることはよく知られており、

偏微分方程式だけでなく境界条件の同定は多くの場合避けられない問題で ある。分布定数系の近似手法や解析手法のほとんどは境界条件に依存する が、 その同定のためには境界条件に依存しない統ーされた近似手法を用い るなどの配慮が必要である13)46)47)。

境界のタイプを例えばノイマン型と仮定したときの未知境界入力を推定 する問題は線形推定問題として取扱えることがある 48)。しかし、システム を入出力関係を表すモデルと見るとき、時間的に変動する未知外部入力を 推定することは実際上困難である。そこで、 境界入力を既知とし境界条件 に現れるパラメータを同定することにより モデリングを行うアプローチが 考えられる。砂原らは2次元領域で定義された放物型分布定数確率システ ムにおいて、入力の存在しない境界における未知パラメータの推定問題を 定常特性に着目して解いたが49)、プラントの過渡状態に関するデータから 推定を行うものではない。一方、 Pierce50)はl次元放物型分布定数システ ムに対し、 可同定性問題を境界のパラメータについても考察しその解を求 めているが、同定手法そのものを具体的に明らかにしているわけではない。

本章では、 プラントの有限時間内の入出力データを用いて境界条件に現 れるパラメータを決定する手法を提案し、そこにおいて実際に必要な計算 手続きを示している。

センサーの数や配置が制限され、観測データ数が限られている状況では 出力誤差法 41)42)による定式化が有利である。 しかし、境界条件は近似モデ ルおよびモデル出力へ複雑に影響するので、この場合同定問題は非線形最 適化問題として取り扱う必要があり、そのためには未知パラメータを適切 に反映する近似モデルの構成が一つの問題となる。

最適化問題においては評価関数がパラメータに対して、微分可能な程度 に滑らかであることが必要であるから、モデル出力がそのような条件を満 たすよう近似モデルを考慮する必要がある。 そこで全パラメータ領域を統 ーされた解析手法で取り扱い、近似モデルの各要素がパラメータに対して 連続であるよう考慮しなければならない。

5. 2 問題設定

l次元領域上で定義された状態u(t，x)，xε[0，1]の時間的空間的挙動が 、 つぎの1次元線形放物型定係数分布定数システムで定まるとする。

au 汀211

云(い)==α

お

^{(い)+ b1} (x)川)， t >仏O<x<l (5.1) ßLU(t，O)一(1-ßL) ð� δ u (t， 0) ==ηLfL(t)， t > 0 (5.2) ßRu(t，l )+(トグR)石(t，1) == 7JRfR(t)， t > 0 θu (5.3)

u(O， x) == uo(x)， 0::; x三l (5.4)

ここで ßL，ßRε[0，1]， α> 0， ηL > O，'7JR > 0とする。(5.1)式は放物型の偏微 分方程式であり、 (5.4)式は初期条件、 (5.2)式および(5.3)式はそれぞれ x == 0とx== 1における境界条件である。境界条件に現れるグL，ßRをそれ ぞれ領域左端、領域右端における境界パラメータと呼んでいる。ßL== 1の ときx==Oにおける境界条件はデイリクレ型、ßL== 0のときノイマン型、

0< グLく1は混合型である。

本章の目的は、この境界パラメータの片方もしくは両方の未知パラメー タ をプラントの観測データから決定することである。そこでプラントから の観測データYda同(ti) に最もよく適合するようなモデル出力似た;q)が得ら れる境界条件のパラメータq== ( グL，ßR)εn == [0，1] ^x [0，1]を、次のような

規範のもとで見いだすよう に定式化する:

出

^{J(q)， J(q):==}

乞

Ily(ti;q)-Ydata(ti) 112

i=l

y(t; q) == [U(t，と1) ... U(tぺm) ]T，

0三c;l.<・・・くとm::; 1。

(5.5)

(5.6) ただしたはt番目の観測時刻である。またqの取り得る範囲Q _{C R}² は_閉 領域であり， これをここではパラメータ領域という。

(5.6)式に示す通り、 モデル出力似たq)は(5.1)-(5.4)式で定まるモデルの状

態u(t，x)の観測位置{ご1，.・・?とm}における点観測によるものとし、初期条件 u。とシステム入力fL，fn， f1 は既知とする。

モデル出力ν(t;q)を解析的に得ることは特殊な場合以外は困難である。そ こで以下では分布系モデルに対してN次元の集中定数系近似モデルを構成 し、これを用いてパラメータについて連続なモデル出力の近似値yN(t;q)，qε

Qを得ることを考える。

5. 3 近似モデルの構成

(5.5)式を最小にするqを勾配法などによる直接的な最小化により実現す るためにはモデル出力ν(t;q)を数値的に求めることが必要である。そこで 本章では， パラメータに依存するモデル出力の近似値yN(t; q) を得るため のN次元の集中系近似モデルの構成を説明する。

重み付き残差法は試行関数の取り方によって内部法， 境界法， 混合法が 知られている。同次境界条件を満足する試行関数を使用する内部法は、 境 界条件を忠実に反映することになる。境界条件の同定問題が偏微分方程式 のパラメータの同定問題に比べると難しいとされるのは境界条件が未知の 場合に取扱いが複雑になりやすいためである。

(5.1)-(5.4)式で表される様なシステムについて3. 2節で与えられる近似 モデルは次のようなものである。

ただし

dU

�; (t) = A(q)U(t) + B(q)f(t) ν(t) = C(q)U(t) + D(q)f(t)

AN := M-1(φ，AφT)，

BN := M-1(φ，A\þT +B) _ANM-1(φ?宙T)，

(;N(X) :=φT(X)，

iJN(X) := \þT(X)ーが(x)M-1(φ，\þT)，

φ(x) := [ゆ1 (x) の(x) ・ ・・ ゆN(X)]T，

宙(x):= [7þ・・・ψm]T M:= (φヲザ).

(5.7)

(5.8)

試行関数ψL(X)，ψR(X)として通常考えられるのは、 ßL= ßR = 0の場合はz

の2次多項式48)、それ以外の場合zの1次多項式に選ぶことである34)0 3.

3節で述べた静特性補償の観点、から一定入力に対するシステムの定常解とし ての1次多項式に取ることが望ましいが、両端がノイマン型(ßL=ゐ= 0 ) の場合には、 そのような1次多項式は存在しない。そこで境界パラメータ

の同定のために全パラメータ領域に対して連続な試行関数として次式のよ うな2次多項式を用いることを提案する。

CPL(X) :== (1-x)2/(2 - ßL)

CPR(X) :== x2/(2 - ßR)

ψL(X)ヲψR(X)の選択の違いにより近似モデルは何らかの影響を受けるが，

3. 2節の命題により、問題の試行関数は座標関数の範囲内の違いでは近 似モデルに影響が及ばないことがわかる。

(5.9)

5. 3. 1 固有関数を用いた近似モデル

正規化された固有関数を用いた場合、 ことなる固有値に属する国有関数 は互いに直交するという性質より行列Mおよび

何

^，AφT

)

^{は次式のように}

非常に簡単になる。

M == 1，

(

^φ，AφT

)

⁼⁼ diag[入1， À2，・ ムN] (5.10)

ただし、 Iは単位行列である。 bL， b R， b Ld， b Rd の各成分はそれぞれ以下 のように計算される。

ßL == ßR == 0の場合

N

1i

一一一一

n4

?N qL

‘，

π

­

qL

・ 114

?

・ 1i nノ“

ト一一一一

ソ i

\・2.2

? /

月々 ? ?

1 1 αo

rlJ、tkrlJ】L

一一一一 R M

' hu rz /fk 'o

一一一一L1川リ

'hu t

/IKL

'O

その他の場合

(bL)i == L3 7t {2α; cos ki -2(1 - a;) sin ki kt(2 -ßL)

+ 2ki(1 - a;) +α;(k; - 2)}

(bR)i == L3 γ�{[ -a;(k; + 2ki -2) + 2ki] cos ki k1( 2 -ßR) + [-a�(k; -2ki -2) + k; -2] sin ki -2α2}

(bLd)t= 2aìi 0 \ {( 1 - an sin ki +α;(1 -cos ki)}

ム(2-ßL)

(bM)t= 2αγt {(1-d)sinkt+α; (1 -cos ki) } ん(2-ßR)

ここで例えば(bL)i ，i == 1，・・・，Nはベクトル bLの第t成分を表す。以上の結 果において 仇→0，内→0の極限は仇==ßR == 0 の場合と一致することが確 かめられる。すなわちパラメータに対するモデルの連続性というパラメー タ同定において要求される一つの条件を満たしていることになる。

5. 3. 2 3次スプライン関数を用いた近似モデル

近似モデルを数値積分によって得るのは全く容易であるが計算コストが 固有関数に比べ著しいため、数式処理を用いて解析的に積分計算を行って おくohj=3，4 ，・・・，N -2についての結果を次式に示す。

151h/315，

397h/1680，

Mi，j ==

�

^{h/4 2，}

h/504 0，

。

2/(3h)，

1/(8h)，

(ゆi，Acþj) ==

�

^1/(5h)，

1/(120h) ，

z-j==Oのとき li -jl == 1のとき li -jl == 2のとき li -jl == 3のとき li -jl三4のとき

z-j==Oのとき li -jl == 1のとき li -jl == 2のとき li -jl == 3のとき li -jl三4のとき

左端に関係する i== 1，2またはj== 1，2の場合については次に示す通りであ

る。右端についてもこれと全く対称である。

Ml=

[;ijj ^d

⁺

会

^LbL+

会

^b

ⁱ ]

^h

I 43 354 1

_

¹

_ _

^I

Ml，2 = M2，1 =

I

商

α

L

C

^L+ ₁^-_6-8� ^{aLdL +}

ず

^L

C

L + 84 b L d L

I

^h

1

=

[会

^L+

命中

M1.4 =川1 =

[お ^αL ]

^h

Mつつ= [土C� "，"

^{I 252 Lt +}

. ^土

^{42 川 +}

LJ LJ • ^坐�

^{8820 d�Lt}

I

^{I h}

同3=同2=

[布 ^C

^L+

品中

崎4=川2 =

[か

^L

]

^h

ぬ5=同2 =

[命中

11一LM

『lilt-」 rU JU

L

'nu

--ω + L

ρしrb 'hU

1一ω+

rb

Jl一hl一hl一h向1||」1llJ

1lJ3一O

k

- -h仇2Ll・-F4

.

LU 一ハU

可l1111J

一nu

1i一 00

- -一つ山1一h2L

.

1 7nl一3一l

d

+ - L

一1lili--一O-1一h1|||」+〆+

仇 2一3 九11J仇 L G

U

一O

一一Oむ

一o

'D 1i一 ハU G 1ょ-nL 1i一つ臼

6

1 4一つ 臼

L l一4 1

70

一ー にも一1170一1

α

一一一一一一ド一副 33-1JL

r il--」rill--Lrill--L G rill--」 「lilt-L FIli--L -

c 一一 一一 一一

.

1

70

+ 0 0 0 + ー が，

d d

2LAVAVAVAYAVAV -fJ

A A A A九 A

AA

1i-っυ ? ?

?

"

lJ

4 1 i - rt 3 4 5

一仇AWAVE AV AV

AV

，，，‘、、 ，，，巴、、 ，，，t、、 FEBBS'''EBEE-L ，，，‘、、 ，，，‘、、 ，，，‘、、

一一一一一一一一一一一一一一一一

14

2 3 4

K 1

3 4 5

AVAVA

V AV -丸 AV AV AV A A A A UW A A A

14

14 V つ剖 内 L

AVAVAYAYL，

q L

φ AV

u 4 fi - -{、、 /{K It-

AV

r

，，{、‘

/{、、 ，，tI ，，，、、、

さらにbL ， b Ld， b R ， b Rd のそれぞれの第t要素は次式の通りである。

(ψL，仇)==

(ψR，仇)==

h /J \{(60αL + bL)h2ー(168αL+ 6bL)h 360(2 -ßL)

+180αL + 15bL}，

h {(CL+479dL)h2ー(6CL+ 726dL)h 360(2 -ßL)

+15cL + 345dL}，

J... 3 ( 1、

ピ

_{L， 一 μL}^ニo

�

_l ^(i-N?+

去

ù J

}

h3 (CR + 479dR)，

360(2 -ßL) h3

(60aR + bR)，

360(2 -ßL)

h (60αL + bL)，

360(2 -ßR)

h (CL+479dL)，

360(2 - ßR)

品 (

^。-1)2+

j)

^?

h {(CR+479dR)h2ー(6CR+ 726dR)h 360(2 -ßR)

+15cR + 345dR}，

h

/J \{(60αR + bR)h2ー(168aR+ 6bR)h 360(2 -ßR)

+180αR + 15bR}，

(AψL，仇)==ァ2h ーア? 15t三N L， 一 μL

(AψR，仇)==一一一2h 1くi< N 2 -ßR'

i==2

3<i<N-2 i==N-l i == N

i == 1 i==2

3<i<N-2

i==N-1 i == N

5. 4 出力誤差法による境界パラメータの同定

3. 2節で構成した 近似モデルを用い解析解(5.27)式から生成した入出力 データにより境界パラメータの同定 を行う。近似モデルの次数は N == 6と した。観測 点の位置はと1 ==♂/6 == 0.28868， Ç2 ==♂/2 == 0.70711 に置く。制 約条件付きでの極小点を 探索を行う のにパウエル法33)を用いた。計算終了 の条件として収束判定値に 10-4、関数評価回数上限に 40000を与えた。評価 関数の導関数はパラメータによるステップ巾10-8の前進差分の差分商で近 似した。ただし、パラメータ領域の境界 上の点で前進差分が取れない場合

(ßL == 1 and/or ßR == 1) ，後退差分に切り替えることとした。

観測時刻をた== i/100， i == 1，・ぺ10としたとき、 初期値を何点かにと って みても、 ßL を正しく推定した ものは なかった (Table.5.1)。これはFig.3.5(a) を見ても分かるように時刻 0.1まででは x==Oの境界 におけるシステムの励 起が十分でないためである。

観測時刻をt == 20/100 まで 2/100ステップ毎に取 ったときは 、N � 6にお いて上の例でのすべての初期値において満足な結果が得られた。初期値を

ßL == ßR == 0.1とした 場合のN == 2，4，6，8 における 結果を Table.5.2 に示す。

また N == 6の場合のモデル出力の計算回数とパラメータの探索値と の関係 をFig.5.1に示す。

Table. 5.1 同定結果(観測時間t==O.lまで)

ßL==l，ßR==l，fL== f1 ==O，!R==sin 9πt ，N==6

グL出発値 内出発値 グL推定値 九推定値

。 。 6 X 10-5 6 X 10-5

。 4 X 10-18 1.0000

。 0.0000 1.0000

1/4 0.2499 1.0000

1/4 3/4 0.2500 1.0000

1/2 0.4997 1.0000

1/2 0.4999 1.0000

1/2 l 0.5000 1.0000

3/4 0.7494 1.0000

3/4 3/4 0.7499 1.0000

1 _。 -1 X 10-16 1.0000

1 0.0000 1.0000

1 0.0000 1.0000

Table.5.2固有関数による同定結果(観測時間t==0.2まで)

真値:グL==ßR==lう初期値:ßL=白=O.l，!L=!l=O，!n=sin 91ft 近似次数N ^{ßL推定値 内推定値} Jmin

2 0.0952 0.9031 6.51 X 10-2 4 1.0000 0.9975 1.30 _X 10-3 6 0.9981 1.0000 1.15 X 10-4 8 1.0000 1.0000 1.63 _X 10-5

Table. 5.3. 3次スプライン関数による同定結果(観測時間t==0.2まで)

グL==l，ßR==l，fL==O，fR==sin 9πt

近似次数N ^{ßL推定値 九推定値} Jmin 6 0.9990 0.9999 3.20 _X 10-6 7 1.0000 0.9999 7.43 X 10-8 8 1.0000 1.0000 1.98 _X 10-7 9 1.0000 1.0000 5.54 X 10-10

ト\勾

担 {同 組0.7

々\

ス0.4

1・ l i --

(3�)

l ' i

l

--

}了\

弓〆 /

ヘ

0.1

0 5 10 15 20 25

モデル出力評価回数t

30

Fig. 5.1モデル出力の計算回数とパラメータの探索値

αコト仏

5. 5 結言

境界入力の存在する線形境界条件を持つ1次元放物型定係数分布定数シ ステムにおける境界パラメータの同定問題を出力誤差法で定式化した。 こ の手法はセンサーの数や配置が制限され、観測データ数が限られている状

況では極めて有効である。 この問題を解くために、未知パラメータを適切 に反映する近似モデルの構成を行い、座標関数として固有関数と3次スプ ライン関数をそれぞれ用いた場合の近似モデルを明確にした。 またこれら の近似モデルを用いて、 システムの励起が十分ならば境界パラメータを同 定し得ることを明らかにした。

第6章 結論

各章の要約をもって結論とする。

第1章では、研究の背景と目的および論文の構成を明らかにした。 第

2章では、本論文で取り扱う放物型分布定数システムを定式化した。また、

重みつき残差法による有限次元モデルの従来の構成について説明し、 それ における問題点をモデリングの観点から指摘した。さらに、分布入力に対 する手法と境界入力に対する手法とは従来別々に考慮されていたが、境界 入力は理論的には分布入力に帰着して取り扱えることを説明した。

第3章では、静特性の解析の結果をベースにして動特性の近似モデルを

構成することを提案した。定常解析の結果を組み込んで近似解を仮定しで も、簡潔な形の近似モデルが得られることを示した。このようにして得ら れた近似モデルの定常特性が、元の分布定数系の近似モデルと確かに一致 することを証明した。さらに固有関数近似においてモデルとシステムのそ れぞれのモード的な解釈を示し、動特性の誤差の性質も理論的に明確な形 で表されることを明らかにした。提案した手法は時間領域でのシミュレー ションだけでなく周波数応答上での誤差の比較でも低い近似項数で優れた 性能をもたらすことを示した。また本手法は分布入力のみならず境界入力 に対しでも、統合的に適用が可能でかつ同様に良好な近似精度を与えてい ることを明らかにした。

第4章では、固有関数を求めることができないような系の近似モデルの 構成を解析的に可能にするための一つの試みを示した。スプライン関数の 微積分演算が数式処理に適していることに注目し、 その援用による解析的 な近似モデルの構成を示す。これによって数値計算の精度や計算負荷の軽 減などが図られ、 パラメータ推定などの種々の問題へ効率よく応用できる ことを述べた。さらに3次スプライン関数が、座標関数として定評のある

固有関数と同等かむしろ優れた性能を持ち得ることをシミュレーションな らびに周波数特性上の比較で定量的に明らかにした。

第5章では、 入力が存在する境界条件に未知パラメータを有する放物型 分布定数系の同定問題を出力誤差法により定式化した。未知パラメータに 対して連続な近似モデルを構成することが必要であり、 これを前章までに 論じた結果を適用して導いた。座標関数として固有関数と3次スプライン 関数をそれぞれ用い 、 システムの励起が十分ならば境界パラメータを同定 し得ることをシミュレーションにより示した。

謝辞

昭和6 2年本学大学院工学研究科に入学し、本研究に着手して以来今日 に至るまで、本学電気工学教室の相良節夫教授には終始懇切なご指導とご

討論を頂いた。 また、本学大学院工学研究科電気工学専攻の長田正教授、

情報工学専攻の島崎箕昭教授には暖かいご助言を頂いた。 ここに謹んで緒 先生方に感謝の意を表します。

また、本研究を行う上でさまざまな貴重なご意見を頂いた本学工学部和 田清助教授、村田助教授ならびに工学部電気工学科制御研究室、工学部一 般電気工学教室、総合理工学研究科エネルギ一変換工学専攻エネルギ一変 換制御工学研究室の先輩諸氏、大学院生の皆様に対しでも心よりお礼申し 上げます。

さらに、本研究において有益なご助言を頂きました福岡工業大学の中野 和司先生、京都工芸繊維大学の大住晃先生をはじめとする皆様に感謝いた します。

参考文献

1)植木， 池田:環境汚染制御[II]，計測と制御， 13-10， 781/796 (1974)

2) ^{W. R.} Grah am: Control Operations in Advanced Aerospace Systems， IEEE

Control Systems Magazine， 7-1， 3/8 (1987)

3)野寺隆:計算を行うためのツール， システム/制御/情報， 35-11， 665/673 (1991 )

4) ^J. M. Boyle， M. P. Ford and J. M. Maciejowski : M山ivariable Toolbox for Use with MATLAB， IEEE Control Syst�ms Magazine， January 1989， 59/65

5)相良，中野:集中と分布[V卜有限次元近似のための数値計算，計測と制御?

26-12， 1065/1073 (1985)

6) G. D. Smith: Numerical solution of Partial differential equations - Finite Diι ference Method， 3rd ed.， Oxford (1985)

7)フオーサイス， ワソー(藤野精一訳) :偏微分方程式の差分法による数値 解法， 吉岡書店(1960)

8)フインレイソン(鷲津他訳) :重みつき残差法と変分原理， 培風館(1974) 9)嘉納秀明: 重みつき残差法による分布系熱交換器の近似，計測自動制御

学会論文集， 12-6， 632/636 (1976)

10)藤川英司:重みつき残差法の適用によるむだ時間系の集中定数系近似，計 測自動制御学会論文集， 15ーム609/615(1979)

11 )市田， 吉本:スプライン関数とその応用， 教育出版(1979) 12)鷲津， 池川:有限要素法， 岩波(1987)

13)相良，中野:分布定数系同定へのアプローチ，計測と制御， 19-11， 1035/1043 (1980)

14) ^{H. T.} Banks : A survey of some problems and recent results for parameter esti­

mation and optimal control in delay and distributed parameter systems， Lecture

Notes in Pure and Applied Math， 81， Eds. K. Hannsgen et al. ， Marcel Dekker，

3/24 (1982)

15) _H. T. Banks， J. M. Crowley and K. Kunisch : Cubic Spline Approximation Tech­

niques for Parameter Estimation in Distributed Systems， 1EEE Trans. Automatic Control， AC-28-7，773/786(1983)

16) _{H. T.} Banks and P. D. Lamm : Estimation of Variable Coefficients in Parabolic Distributed Systems， 1EEE Trans. Automatic Control， AC^-30-4_，386/₃98(1985) 17) _{H. T.} Banks : Spline Approximation for Functional Differential Equations， Jour­

nal of Differential Equations， 34， 496/522 (1979)

18) H. T. Banks and K. A. Murphy : Estimation of Coefficients and Boundary Param­

eters in Hyperbolic Systems， SIAM J. Control and Optimization， 24-5， 926/950 (1986)

19)斉藤制海:数式処理雑感(2)-数式処理の工学分野への応用?システム/情 報/制御， _33， 600/601 (1989)

20)伊藤清三:偏微分方程式， 培風館(1966) 21)伊藤清三:拡散方程式， 紀伊園屋(1979)

22) H. O. Fattorini: Boundary Control Systems， S1AM J. Control， 6-引1968) 23) _{A. V.} Ba叫la比kr巾i治shr凶n: Applied F\町lnc

24判4め) ^{R. F.} Curtain and Pritchard : 1r凶n凶凶it旬e Dimensional Li凶nea紅.r Systems T、heoαr、

Springer (1978)

25) J. S. Gibson and 1. G. Rosen: Approximation ofDiscrete-time LQG Compensators for Distributed Systems with Boundary Input and Unbounded Measurement， Au­

tomatica， 24-4， 517/529 (1988)

26) _{J. R.} Cannon， S. P. Esteva加d J. Hoek: A Galerkin Procedure for the Diffusion Equation subject to the Specification of Mass， S1AM J. Numer. Anal.， 24-3，

499/515 (1987)

27) 1. Lasiecka : Ga叫ler比ki∞i

Problems Wi比th Rough Boundary Data一Lらp Theory， Mathematics of compu­

tation 47-175， 55/75 (1986) 28) R. Vichnevetsky: U se of Func

Soluもtion of Initial Value Parもtia叫1 Differenもtia心1 Equation Problems， IEEEτ'rans.

Computers， C-18-6，499/512(1969)

29) T. 1. Seidman: Approximation Methods for Distributed Parameter Systems，

in Distributed Parameter Control Systems， Theoryαnd application， Ed. S. G.

Tzafestas(1982)

30) M. R. Cl刈a山ara and E. J. Davison: Correspondence on a Method for Simpli­

fying Liear Dynamic Systems ; IEEE Trans. Autom. Control， AC-12， 119/121 (1967)

31)国松， 浜田:集中・分布システムの安定論， 実教出版(1988)

32)嘉納秀明:分布定数系の集中系近似モデル，計測と制御， 19-11， 1044/1050 (1980)

33)富士通: SSLII使用手引書(1987)

34)ファーロウ(伊理訳) :偏微分方程式; pp.31-51啓学出版(1982)

35)相良，今井:境界入力を持つ放物型分布定数系のガラーキン法による近 似モデルについて， 計測自動制御学会論文集，27-12， 1366/1373 (1991) 36)今井，相良:放物型分布定数系のガラーキン近似法における静特性の補

償について， システム制御情報学会論文誌(投稿中)

37) K. Hara， T. Yamamoto and K. Ougita : Parameter identification of distributed parameter systems using spline functions， International J ournal of Systems' Sci­

ence， 19-1， 49/64 (1988)

38) P. M. Prenter : Splines and Variational Methods， Wiley(1975)

39) H. T. Banks and 1. G. Rosen : Fully Discrete Approximation Methods forもhe Estimation of Parabolic Systems and Boundary Parameters， Acta applicandae Mathematicaβヲ7， 1/34 (1986)

40)大野義夫:U，niform Bスプライン曲線の重み関数， REDUCEプログラミン グ資料第三集， p.413， 東京大学大型計算機センター(1986)

41)古田勝久:線形システムの観測と同定;コロナ社， p.158 (1976)

42)相良， 秋月， 中溝， 片山:システム同定;計測自動制御学会， p.8(1979) 43)相良節夫: モデリングとシステム同定， 電気学会論文誌C分冊， 108ふ

292/297 (1988)

44) M. P. Polis and R. E. Goodson : Parameter Identification in Distributed Systems:

A Synthesizing Overview; Proc. IEEE， Vo1.64-1， 45/61 (1976)

45) C. S. K山rusly : Distrib凶ed param伽system identification A survey， Int. J Control， 26-4， 509/535 (1977)

46) S. Aihara and A. Bagchi : Parameter Identification for Stochastic Diffusion Equa­

tions with Unknown Boundary Conditions， Applied Mathematics and Optimiza­

tion， 17， 15/36 (1988)

47) K. Ohnaka and K. Uosakai : Boundary element approach for identification of boundary conditions of distributed-parameter systems， Int. J. Control， 41-4，

981/990 (198S)

48) _H. Miyauchi， S. Saga肌M. Maeda and J. Imai : Estimation of U nknown N on­

linear Boundary Condition in a Parabolic Distributed Parameter System; Proc.

IMACS/IFAC

Symp. Modelling & Simulation of Distributed Parameter Systems， 73/78 (1987)

49) _Y. Sunahara， F. Kojima et al.: Identification of Boundary Conditions for a Two Dimensional Diffusion System under Noisy observation， Prepr. 17th JAACE SSS 205/208 (1985)

50) A. Pierce: Unique Ide凶五cationおEigenvalues and Coefficient in a Parabolic Problem ;SIAM J. control & Optimization ，Vol.17-4 494/499 (1979)

付録 ^A. 1次元系に対する固有関数近似

ガラーキン?去における座標関数仇(めにはシステムの固有関数がよく用 いられる。システム (2.6)-(2.8)式の固有関数の一般的表現は次のようになる

ことが知られている。

仇(X):==γi(αi2 sin kix ⁺ (1 - Ct:i2) cos kix)

ん- αki2 ， i == 1， 2，・・・，N (A.1)

境界パラメータグL，ßR についての固有関数を得るには， (A.1)式が 次 式の同 次境界条件を満足す るように αt2?んを定めればよい。

ßLØi(O)一(1- ßL)

学

^例==0

12 (A.2)

ßRØi(1) ⁺ (1 - ßR)

去

^(1)==。

以下ではその結果を示す。

i) ßL == ßR == 0のとき

(札仇) == 1， i == 1ぃ・・，N

わ:==(i-1)久 向:==0， i==1，.・・，N

i == 1

i==2ぃ・・，N ii) ßL == ßR == 1のとき

(A.3)

(A.4)

ん:==zπ? αi :== 1， γf:=2? t=lγ・・，N (A.5) iii)その他の場合

ん:区間((i- 1)π?付)上の(5.23)式の一意解 注)ここでの勺は2. 2. 2節に現れるγとは無関係である。

{ßLßR

-(1 -

ßL)(l

-ßR)k九sin ki

+

{ßR(l -

ßL) + ßL(l

- ßR)}ki cos ki

=

0 (A.6)

ぷ- ^/\

FL

ⁿ )ki^' (A.7)

f =(l+l ーグ

2 ' 4kィ sin 2ki

α汽l一 向2) /___O\L ^{， 0\ 1} '1\1 -1

- 2k

t (ω2ki ⁺ 2ki - 1)

J

^(A.8)

i)，ii) の場合については明らかである。iii) についても(A.1) 式を(A.2) 式 へ代 入して整理することで(A.6)，( A. 7)式が容易に得られ， (A.8)式も正規化の条 件(A.3)式 から直接求まる。 モードの決定方程式(A.6)式 についての閉区間 ((i - 1)π?付)上での解の存在と一意性(以下を参照)により、 その数値的な 求解は挟みうち法などを用いて容易に行える。

決定方程式(A.6)の解の存在範囲

解の存在範囲とその一意性について説明する。j を負でない 整数として kj i= j7rとする。(A.6)式の両辺をsinんで除くと次を得る。

tan -1 k � ^._，

= _

^. --.-_(_1 - ßL)(l - ßR) ι

ßR(l -

ßL) +グL(l

- ßR) 叶

FLßR

¹

ßR(l -

ßL) +グL(l

- ßR)ん

左辺はんの関数として区間((i - 1)7r， 付))で ∞から一∞への単調減少であ り右辺はα> 0， ^b > 0 をあるスカラとしてα丸一αjkiの形をしているから全 区間に渡って単調増大である。従って， (5.23)式の解は区間((i - 1)π?付)で一 意に定まる(Fig.A.1 ) 。

段(も税制時択限定(もり{ll卯'

一〈・切一』

付録 B. 3次スプライン関数近似におけるREDUCEソースリスト

本ソースはREDUCE3.3により、近似モデルの積分行列の計算結果をFOR­

TRAN形式で出力する。 リスト中に現れる変数Pは本文中のM ==

(

^φ，φT

)

に、変数Qは

(

^<Þ，AφT

)

^{にそれぞれ対応する。}

‘�'l.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.'l.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.'l.'l.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.'l.Y.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.

Y.Y. スプライン近似計算プログラム

Y.Y. ( 1次元放物型分布定数システム版)

'l.Y. 1 9 9 0年2月完成

Y.Y. 作成者 ^: 中城康平 ^・ 今井純

、�Y.‘

‘�Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.'l.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.Y.

on fort$ Y.FORTRAN出力 on gcd$ Y.既約分数

out "e:\result"$ Y.出力ファイル名

‘�Y.Y.Y.Y.シス子i入漁填=干の定義、�Y.Y.'l.Y.Y.Y.Y.Y.Y.Y.Y.Y.

procedure system(u);

df(u，x，2)$ Y.円柱座標Z方向J (定係数系) Y. df(u，x，2)+1/x*df(u，x)$ Y.円柱座標半径方向

Y.Y.'l.Y.Y. obh = 1/(6申h彬申3) Y.Y.'l.Y.Y.Y.Y.Y.Y.Y.'l.Y.

Y.h := hnode$

Y.Y.Y.'l.Y.スブラベγ現数の支義‘�'l.Y.'l.Y.Y.Y.Y.Y.Y.Y.'l.‘

y.y.y.y.y. bb(1，jj):区間[x_{i-3+jj}，x_{i-2+jj})におけるB-i(x) y.y.y.y.y.y.y.y.y.y.y.y.y.y.y.

matrix b(1，4)$

matrix bb(1，4)$

b(1，1) := obh事(x-xm2)ホ申3$

b(1，2) := obh*(h料3 + 3本h**2市(x-xm1) + 3*h*(x-xm1)材2 - 3事(x-xm1)料3)$

b(1，3) := obh申(h料3 + 3申h*申2申(xp1 - x) + 3*h事(xp1-x)村2 - 3*(xp1-x)料3)$

b(1，4) := obh事(xp2-x)料3$

for jj:=1 step 1 until 4 do<<

bb(1，jj)

:= sub(xm2=(ii-2)内+ll，xm1=(ii-1)*h+ll，xp1=(ii+1)内+11，xp2=(ii+2)内+ll，b(1，jj))

>>$

y.y.y.y.y. ^aL bL cL dL aR bR cR dR (境界上の基底結合係数 )の計算 y.y.y.y.y.y.y.y.y.y.

matrix a(2，2)$

matrix work(2，1)$

'W''W' := 0$

indefl := int(sub(ii=0，bb(1，3))，x)$

'W''W' := 'W''W' + sub(x=l吋l・t11，indefl) - sub(x=ll+O，indefl)$

indef2 := int}sub(ii=0，bb(1，4))，x)$

'W''W' := 'W''W' + sub(x=2申h+11，indef2) - sub(x=1*h+11，indef2)$

a(2，1) := 'W''W'$

'W''W' := 0$

indefl := int(sub(ii=-1，bb(1，4))，x)$

'W''W' := 'W''W' + sub(x=l吋t+11，indefl) - sub(x=O+ll，indefl)$

a(2，2) := 'W''W'$

a(1，1):=beta1申sub(x=0+11，ii=0，bb(1，2))ー(1-beta1)事sub(x=0+11，ii=0，df(bb(1，2)，x))$

a(1，2):=beta1*sub(x=0+11，ii=-1，bb(1，3))ー(1-betal)*sub(x=0+11，ii=-1，df(bb(1，3)，x))$

Z 'W'ork

:=l/a場mat ( (0) ， (h) ) $ x

a1:=sub(obh=1/(6申M車3)，'W'ork(1，1));

b1:=sub(obh=1/(6事h料3)，'W'ork(2，1));

'W''W' := 0$

indef := int(sub(ii=-1，bb(1，4))，x)$

w'W' := 'W'w + sub(x=l吋l+ll，indef) - sub(x=O+ll，indef)$

a(2，1) := w'W'$

'W'w := 0$

indef := int(sub(ii=1，bb(1，2))，x)$

'W''W' := 'W''W' + sub(x=lぬ+ll，indef) - sub(x=O+ll，indef)$

indef := int(sub(ii=1，bb(1，3))，x)$

'W''W' := 'W'w + sub(x=2吋t+11，indef) - sub(x=l吋t+11，indef)$

indef := int(sub(ii=1，bb(1，4))，x)$

w'W' := W'W' + sub(x=3吋t+1l，indef) - sub(x=2吋t+11， indef) $ a (2 ， 2) : = 'W'w $

a(l，l) :=beta1句ub(x=O+ll，ii=ー1，bb(1，3))ー(l-betal)句ub(x=0+11，ii=-1，df(bb(1，3)，x))$

a(1，2):=betal*sub(x=0+11，ii=1，bb(1，2))ー(l-betal)本sub(x=0+11，ii=1，df(bb(1，2)，x))$

1 'W'ork

:=l/a本mat((O)，(h))$

X

c1:=sub(obh=1/(6*h料3)，'W'ork(1，1));

d1:=sub(obh=1/(6事h料3)，work(2，1));

、‘，，，、‘，J、、.J、‘.J

14 1ム 1占 1ム ab c d

uuEU

+LV +lu +b +b e e e e bbbb

=

=

=

= 1 1 1 14 a a a

a

+LW +lu +lw みlu e e e e

bbbb ，f、，，E、，，E、，，E、bbbb

u u u u qu 冒M EM -U z --

=

= r r r r

ab

e d

'l.'l.'l.'l. a1，bl，c1，d1，ar，br，cr，drを境界パラメータbeta1，betar 'l.'l.'l.'l.'l.'l.'l.

'l.'l.'l.'l.で展開するときは、次の行をコメントアウトする 'l.'l.'l.'l.'l.'l.'l.

c1ear a1，b1，c1，d1，ar，br，cr，dr;日境界パラメータで展開しない

‘'I.'l.'l.'l.'l.'l.'l. 京事与のBス7ラベγ執致 、'I.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.

‘'I.'l.'l.'l.'l.'l.'l.‘ phi1(s，t): IRI可MH(t-l)*h+aλ尊い砧与のキ_8‘'1. 'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.'l.