論 文
Chebyshev級数によるWangerin N(v)
微分方程式に関する二点境界値問題の数値解
相川孝作
高原幹夫 津金修三 (昭和52年8月12日受理)Numerical Solution of the Two Point Boundary Value
Problem for Wangerin N (v) Differential Equation
in Chebyshev Series
KosakuAIKAWA MikioTAKAHARA ShuzoTSUGANE
Abstract The axisymmetric potential problem for a charged disc or sphere in a concentric hole of a gro皿ded sheet is solved by use of Wangerin functions in flat・ring cyclide coordinates (u,v,ψ). The Lapalace equation is R・separated into two Wangerin differential equations, M(のequation and 2V(v)equation. The potential function and the capacitance between two electrodes are expressed by Wangerin functions. The IV(のfunctions, however, do not converge at an end of the interval in this problem and the numerical values can not be obtained. Then we try to get the bounded solution of」N(のWangerin differential equation in a finite Chebyshev series and we can obtain the numerical values of the potential function and the capacitance. 1. まえがき 且at・ring cyclide座標系(t,v,ψ)においてLa− place方程式をR分離するとき,二つのwangerin微 分方程式が得られる1)2)3)4)。その一つの〃に関するN (の方程式は特異点のまわりに級数解として二つの特 解Sp。r(のとTPnくv)を持ち,その一般解は二つの特 解の一次結合としてN(の一CISp。’(v)十C2Tp。’(ので 与えられる。 無限平面導体の孔に置かれた円板または球導体の静 電容量を求めるとき,領域の一端v−Ktで, Sp。t(の, Tp。’(v)が発散して数値解が得られないことが起こ る5)6)。この境界値問題において,N(のを有限個の chebyshev級数の和として表すことにより, N(v)の 数値解を求めることができ,静電容量の数値を求める ことができる。この論文ではchebyshev級数を静電 界境界値問題へ導入する過程について精細に述べる。 2.Flat・ring cyclide座標系とWangerin微分方程式 *日本電気株式会社伝送通信事業部 **昭47年度11月電気学会電磁界理論研究会 Flat・ring cyclide座標系は図一1に示す座標系を縦 軸のまわりに回転したものである。Laplace方程式 は,軸対称の場合には 7・φ一∂㌶d。v{dn〃☆(÷・一怨) +・n・☆(dnv∂φA ∂v)}一・・ O≦u:{9K,0≦v≦K’ (1) で与えられる。ここにαは焦点距離であり,また
一90一
図一1 Kノ 22=(1−sn2udn2v)(dn2v一ん2sn2μ), A=1−dn2usn2v である。なお,snuはsn(u, le)の略, snvはsn(v, ん’)の略(le, ktはそれぞれ楕円関数の母数と補母数) であり,dnu, dnv…についても同様である。またK, K’はそれぞれ完全楕円積分,補完全楕円積分である。 (1)の解を φ=ノdi/2M(u)N(v) (2) とおいて,R分離すると次の二つの微分方程式が得ら れる。 袈+c雲”・誓+(le2・n・u一α)M−・,(・) 際一ゐ’漂〃・漂+(一・n・v+α)N−O,(・) ここにαは分離定数である。 式(3)においてsn2u−z,式(4)においてdn2v−zとい う変数変換を行うと 蔑+丁(⊥+2+2zz−1z−c)筈 2z一βc Z−0 (5) 十 4z(z−1)(z−c) という共通の式が得られる。ただし,式(3)に対しては c−1/k2,β一α,式(4)に対してはc一ゐ2,β一α/le2で ある。式⑤は0,1,cおよび。。の四つの確定特異点を もち,これら確定特異点近傍の級数解がある。式(3)∼ (5)はWangerin微分方程式といわれるが,式(3),(4)を Jacobi形,式(5)を代数形と呼ぶことにする。 いま級数解をu,vの関数として表し,元=0,1,2… とすれぽ,式(3)の特解は S。(le, z)=・= S。(k, sn2u)一Σcゴsn2ゴ・u, (6) 」=O T。(k,2)−T。(le, Sn2U)−S。(le, Sn2めln(Sn2U) +Σ4ゴsn2∫η, (7) ゴニ ここにcゴは展開係数であり,またd」はcゴを決定方 程式の根rで微分してr−0とおいたものである。 d」一芸L。 式(4)の特解は Sρ’(1/le, z)=Spt(1/ん, dn2v) oo =(dn2v 一 1)1/2Σら(dn2v 一 1)グ,(8) ゴ=O Tp’(1/le, z)=7「p’(1/le, dn2v) oo =Σら(dn2v−1)ゴ (9) ノ=0 で与えられる。式(6),(2),(8),(9)はWangerin関数 といわれる。本論文ではこれらをS.(u),Tp(za), Sp’(v),τ〆(ので略記することにする。 snu・・=・Oの場 合,Tp (u)はln(0)=一〇。を含むので,これを捨て る。したがって式(3)の解はSp(n)となり,式(4)の解 は式(8),(9)の一次結合 2v(の一c、S。’(の+c,T。t(の ao) で与えられる。なお,Sp(u), Sp’(の, Tpa’(のをそれ ぞれ図一2,3,4の破線で示す。図の実線は固有関数で ある。 3 三 ↑2 1 0 一1 λ。−0.3536147 / ’ λ1 3.979735 ノ214.02561 ,/ ノ I ’ m
篶霊
d4//__4Zイ1/ 0 λ1 ∼α=0 、、 、 κ/2 κ λ4 一→μ κ2=α5 2 ε亨1
↑ 0 図一2 λ0 1.457000 λ1 7.199171 λ2 18.68349 λ, 35.90997 λ, 58.87861 、ん。//
図一3昭和52年12月 山梨大学工学部研究報告 第28号 3
三2
心1、
0 一1 λ0 0.7392316 λ1 3.615677 λ2 12.22358 λ3 26.57896 λ4 46.67652/
λ1ミこ・−o/
/
/ λ。 →砂jγ2 λ4 図一4 4. 電位関数と静電容量 内側電極をv−Vl,外側電極をv−v2とし,電位を v=Vl,φ=v v−v2,φ=0 とする。この電位関数φは oo φ=(1−dn2usn2v)i/2 Z・AnSPn(u) n=0 {s・・’(v)一:i潟Tw(の} ⑪ ⑫ An= {s・・t(Vl)一享鵠}τ〃ω} ・∫1{SPn(u)}2snudu ⑭ 次にε ==誘電率と置くとき,基準化された静電容量 c/aは c 4πε κ snudnvi dn2usnvicnvidnvi で与えられる。ここtlC SPn(u)は式(3)に境界条件 u=O,M(u)=1, M’(u)=0 ⑬ u=K,M(u)有界,ルf’(の一〇 を与えたときの固有関数である。式(4)のαに式(3)の固 有値を代入したときの式(4)の級数解がSp。’(の,TPn《の である。したがってSp。くのとTp。・(v)は固有関数で ない。よって図一3,図一4の破線で示すようにv−K’で 発散する。またA。は次式で与えられる。 v∫1(Sp。(u)snu1−dn2usn2Vi)du a。∫。(1−、n,u、n,v、)[1−d。,u、㎡Vl ・急繊(u){s・・,(の一鶏器丁〆(Vl)} 一曇。繊ω{s・・”(Vl)一芸{霧・T…’(v・)}]du ⑮
で与えられる。 電位関係,静電容量および固有値の求め方について は文献(5),(6)を参照されたい。 さて,平面導体の孔に置かれた円板の場合はVl=O, v2→K’,同じく球の場合はVl−K’/2, v2−K’である が,式a2)あるいは式⑭に含まれるSヵ。’(〃2), Tp。’(v2) はいずれも固有関数ではないから,v2−K’で発散す るので数値解が得られない。 そこでSp。’(v), Tp。,(v)のかわりにz=ゐ2(v=K’) で展開した関数 Sヵ。’,k・(v), Tρ。’, k・(のを用いてみる と,円板の場合にはv−0で関数が発散するので,結 果としては同じである。球の場合にはle>0.5の範囲 内ではv−K’/2で関数が収束するので数値計算が可 能となる。しかしleの全域については求められな い5}6)o この場合の電位関数は次式で与えられる。 ooφ一Al/2Σ&sヵ。(u)SPn∼た・(の カ=05.Chebyshev近似の適用
⑯ 5, 1Chebyshev近似の概要 級数解が収束しないとき,これを有限個のcheby・ shev級数で表すことにより有界な解を求めることが できる。式a2)のSp。くV), Tp。’(V)を個個に計算する代 わりに,電位φを φ一A’ノ2 Z A。SPn(め」Vp.r(の ⑰ n=O とおいてN(のそのものを直接chebyshev級数によ って表そうと試みる。境界条件はN(0)またはN(K’ /2)が有界,N(K’)が0である。これに占部氏の方 法7)8)の二点境界値問題の解法を適用する。まず,占 部氏の方法の概要を述べる。 領域[−1,1コで有限個からなるchebyshev級数 アm(’)=Σ e。a。T。(の, n・・O,1,2,…〃z ㈹ n=O を考える。a。, al, a2,…a。は未知の係数であり, e。 については皐i三il鷲
である。またTn(t)はchebyshev多項式で Tn(COSθ)=COS nθ (19) である。式⑱からanは an一サ∫lf・(…θ)…nθdθ で与えられる。 fm(t)が領域[−1,1]でtについて連続関数であ るとして,微係数からなる有限個のchebyshev級数 m−1 アm’(め一Σ e。a。’T。(の, n=0,1,2…,〃z−1 n=O ⑳一92一
を考える。ただし a・’一=− P}∫lf・’(…θ)・…θdθ chebyshev級数の直交性とParsevalの等式から式 ⑱のanと式⑳のa。’の問に次の関係が得られる。 an’一Σe、.。9、.。sα、, n−0,1,2,… ⑳* n;O g、−n == 1−(−1)s−n・==O (s−nが偶数のとき) −2 (s−nが奇数のとき) 5.2 平面導体の孔の中の円板 式(4)に対して基準化した境界条件 N(0)=1,2V(Kt)=0 ㈲ を用いる。式(4)は次の二つの連立方程式 N「(v) ・=Y(の 倒 le’2snvcnv Y(v)十(dn2ρ一α)N(の 24 Y「(の一 dnv に置きかえられる。vの領域は[0, K’]であるから次 のような変数変換を行う。
t一芸一1 囲
ただし 一1≦t≦1。 式23),倒はτに関する次の式に置きかえられる。 N「(の=Y(t) 陶 Y’(t)=R(t)Y(t) 十Q(t)N(の 27) ここに R(の一.亙”’2sn{Kt(t+1)/2}・n{K「(t+1)/2} 2 dn{K’(t十1)/2} L灘㌻㌧燃、ヨ ⑳ Q(t)==(与)2{d・・ K’(i+1)一α} ⑳ である。境界条件,式⑳は N(−1)一α1−1,N(1)=・ev2・=O Bo) となる。ここでN(の,Y(のを有限個のChebyshev級数
Nm(の,γm(ので近似し次のように置く。 2Vm(の一Σ e。a。Tn(t), ⑪ n =O Ym(t)一Σ e。b。T。(の β2) n=O ここに,a。, bnはそれぞれm+1個の未決定の係数で ある。 式㈲,B2)をtについて微分して m−1 Nm’(の一Σ e。an’T。(の, B3) n=O m−1 γ仇’(の一Σ e。b。’T。(の. ㈱ n=O a。tC b。’はそれぞれm個の未知の係数である。合計, 未知数が2(m+1)+2m個ある。 式⑳,閻,Bo)をChebyshev近似で表すと 2Vm’(t)==Ym(の, Ym’(t)=R(t)Ym(t)十Q(t)2V.(t), 1Vm(−1)=1,2V肌(の=0 になる。 R(t),Q(t)をChebyshev級数に展開し R(の一Σ e。R。Tn(’), n=O oo Q(’)一Σen(?。T。(の π=0 とおき, のようにChebyshev級数で表す。 R(のγ皿(の一 Z e。rnTn(の, n;O Q(t)AV.(t)−Z]. e。q。Tn(の. n=O 倒 陶 ㈱ BS) Bg) さらに式鯛のR(t)Ym(の, Q(t)Nm(t)を次 ㈹ (4D 式㈹,B2)に対する決定方程式は次のようにして導か れる。すなわち式田,B2), B3)から ant=bn, n=0,1,2,…, m−1, 式㈹,㈱,㈹,㈹から bnLrn十qn, n=0,1,2,…, m−1, 式助から Nm(−1)一Σ(−1)ne。an一α1−1, tl=・O m Nm(1)一Σ e。an一α2ニO n=O が導かれる。ただし a。LΣθ、.。9、.。sas, s=O m b。’一Σe、.。9、.。sb、, s==O n=0,1,2,…,m−1 閏 ⑬ ⑭ ㈲ ⑳「 ㈹ であり,qn, rnは次のようにして求められる。すな わち式⑪において oo m Q(’)AV.(t)一Σ e。Q。T。(の]Z] e、a、T、(の r=O s=O 一丹吋{T・+・(の+T・・一・)(t)}Q・as 一去蜘砺+些。e・2Q・as}T・(t) +婿1自。(en−・e・Qn−・+en・se、Qn・、 *文献(5)の臼の右辺はΣes.。9s.nSasに訂正,またその下の n=o fs,nはes−nに訂正。 { +es一頑・一泌}Tn(t)昭和52年12月 山梨大学工学部研究報告 第28号 であるから qn一彊。(en−・e・Qn−・・+en・…Qn・s+es−ne・Q・−n 十e−ne−sQs)as, 働 n=O,1,2,…,m−1 同様に ㌍鍵。(en−…Rn−・+en+…Rn・・+e・−ne・R・−n 十e−ne−sRs)bs. ㈲ さて2(m+1)+2m個の未知数an, bn, an’, b。’ に対して,決定方程式幽∼㈲において方程式の数は式 92),⑬でそれぞれm個,式⑭,㈲で2個,式⑳’,㈹ で2m個,合計2(m+1)+2m個あるので,これらの 係数を求めることができる。bnが求められれet“ a.が 求められ,式B2),⑪からそれぞれYm(の, AV.(のが 決定される。 まずanとb。との関係式を求めよう。式⑳,⑭から anをbnで表すと an−☆(bn−一d・bn+1)・n−・,2・・…, m a・) ここに
叶㍑::二1隠㌫
である。式㈲を式㈲に代入してa・・−2+b・一†b・+2貢吉器万ぴ 6◎
を得る。したがって係数a。(n−0,1,2,…,n)はbn を用いて式kg),60)で表される。式陶において左辺 y石(t)の係数 bn’は式q6)の関係からbnで,右辺 、R(t)Y(t)の係数rnは式幽によりRn, bnで表され る。Q(t)N(のの係tw q。は助によりQ。, a。で表さ れるが,式kg),60)の関係があるので,結局Q。, b。で 表される。すなわち q。・・一・Qn+丁(L+t(en−iQn−i+Qn・・ { +e・−nQ・−n)}b・ { +e・−nQ・−n)}bi+−
狽p.+9(en−・Qn−・+Qn・・ +墓(、≒}1≒1)Qn+4(、11) ×(en−s−iQn−s−i十Qn+s+1 1 十Us+1−nQs+1..n)− 4(s−1) (en−s+iQn−s+1 +Q・+・一・+e・一・−nQ・−i−n)}bs n==O,1,2,…,m−1 60 となる。 これらの関係式を式B6)に代入して, bs(s=O,1,…, m)に関するm個の一次連立方程式m
Σ、A。、bs−.B。, n−0,1,…, m−1 5=0 を得る。ここに An・一†(瓦+払+÷(en−1Qn−1+Qn・1 十el−nQi−n), 1 1 4 8 +†(en−・Rn−1+Rn・1+el−nR・−n) −es+1_nQs+1_n, (−1)s−1 1 61 A。1−一一Qn+一(e。−2Q。−2+Qn+2+e2−。Q2−。) A・・ =‘(s−1)(s+戸Q・+4(s+1)(en−・一・Qn−・−i 1 十Qn+s−1十es+1−nQs+1−n)− 4(s−1) X(en−s+iQn−s+1十Qn+s−1十es−i−nQs−1−n) +†(en−,Rn−,+Rn・・+u、−nRs−n) −es_ngs_nS s=2, 3,… , m−1, A。m−†(R。・。+Rm−n)−9m−nM, Bn=−Qn, n=0,1,2,…,m−1. 境界条件の一つから得られる式㈲のanに式醐,⑳ を代入して,bsについて整理するとm
Σ Am、’bs−Bm’, 63) s=O ここに AmO’ = 1,、4.i「=1, A・・’一一( γ8二1s−1)(s+1)(s−2・ 3・・… m−1) Ammt=0, Bm’一一α1+α2−−1+0−−1 となる。 式52),63)はm+1個のb、(s−0,1,…,m)に関す る(m+1)個の連立方程式である。bnが決まると Nm(のの係数anは式働,60) rcより得られる。決定方 程式幽∼㈲は式62),63),㈲および60) ae帰せられるわけ である。bnからN.”(t), anから1Vm(の,したがって N’(t),N(のが決定され,式㈱の変数変換を用いて 」V(の=N(の, ⑭ N’(v)一÷N’(の一÷γ(t) 6s) を得る。N(v)が決まれぽ式⑰から電位関数φが求め られ,またNt(のが決まれば基準化容量c/aが求め られる。電位関数は次のようになる。一94一
1.0 ε0.5 ミ ↑ 1.0 ミ ↑… λ2 λ3 K2=0.5 Kノ=1.8541 Kt/2 ■一一一ウv 図一5 K’/2 図一6 K2=0.5 Kノ==1.8541 1,0 ε 2・o・s τ 1.0 ミ ↑… 股72 λ1 λ2 K2=0.5 κノ=1.8541 Kノ/4 −一一一一xv 図一7 K’/4 図一8 K2=0.5 Ki=1.8541 φ一A’/2 Z A。Sp。(u)Np。’(の, n=0 ここに An一γ∫:s九(u)・nudu/凡M・)・ 5θ 6カ 2V(の,φを図示すれぽそれぞれ図一5,図一6のよう になる。静電容量については文献(5)を参照されたい。 5.3 平面導体の孔の中の球 円板の場合と異なるのは境界条件と変数変換であ る。境界条件を 2V(Kt/2)=1, N(K’)=0, 68 変数変換を 4v t=『云7−−3・ ただし 一1≦τ≦1 69) とすれぽ,5.2と同様に1Vρ。’(v)が求められる。電位 関数は co φ一Ai/2ΣA.nSPn(u)Np。’(の, n=0 ここに A・・一・V(1+le)党講鵠埠d・/ N・N・・’(K’2) ㈹ である。式㈹は式06)のSp。’, k・(のの代わりにNp。’(の とおいたものに外ならない。Np。,(v)とφをそれぞれ 図一7と図一8に示す。静電容量については,5.2同様, 文献(5)を参照されたい。なお,計算においてcheby・ shev近以の項数は51とした。 6. 考 察 球の場合のφは,z=カ2(v・−K’)で展開したSPn,, k2 (のを用いた式㈹,またはChebyshev近似から得ら れたNp。くのを用いた式㈹により得られる。 Sp。’, k・ (のとNp。’(のの数値を倍長計算により比較すると有 効数字12桁まで一致する。 また式⑳のR(t)とChebyshev近似による式㈲の Rm(t)の数値計算結果を比較すると,例えば R(t)−0.139939698463874に対して IR(t)−Rm(t)1=0。6613633×10−16 程度の差異に過ぎない,よって計算結果はかなり信頼 できるものと思われる。 謝 辞 Chebyshev近似の微分方程式への応用は故九州大 学占部教授により提起されたもので,本学栗原光信助 教授から紹介され,極めて懇切な御指導をいただいた。 ここに深甚な謝意を表する。 文 献 1)P.Moon and D. E. Spencer:Field theory hand− book. Springer Verl, Berlin(1961). 2)相川,菅沼:Wangerin微分方程式の固有解信学論 55−A,10,pp.564−565(昭47−10). 3)K.Aikawa and M. Takahara:Wangerin functions, Rpt. Fac. Eng. Yamanashi Univ. No.26(1975). 4)相川,高原,久本:Heine関数, Wangerin関数と
昭和52年12月 山梨大学工学部研究報告 第28号 その固有関数の数値計算,電学論A:96,6,pp.287− 294(昭51−6). 5) K.Aikawa and S. Tsugane:The capacity of a charged disc or sphere in a concentric hole of agrounded sheet, Letters in applied and engi・ neering sciences.4pp.257−268(1976). 6)相川,津金:平面導体の孔に置かれた円板または球導 体の静電容量,電気学会電磁界理論研究会資科EMT− 74−11(昭49−8). 7)Urabe, M. Numerical solution of multi・point 8) boundary value problems in Chebyshev series− Theory of the method, Numer. Math. g, pp.341 −366 (1967). Urabe, M.:Numerical solution of boundary value problems in Chebyshev series・A method of computation and error estimation, Lec. Notes Math.109, Confernce on the numerical solution of differential equations held in Dundee, Scot・ land, June,1969, Berlin, Springer Verl. PP.40− 86 (1969).