第7回環瀬戸内応用数理研究部会シンポジウム2003

全文

(1)第 7 回環瀬戸内応用数理研究部会シンポジウム 2003 プログラム・目次（各講演時間には質疑応答を含みます．） 1 日目午後（5 月 6 日火曜日） 15:00. 開会の挨拶仁木滉（岡山理科大学総合情報学部情報科学科）実行委員長挨拶中島正治（鹿児島大学理学部数理情報科学科）. 第１セッション（15:15-16:55） 15:15-15:40 I+R 型前処理付 Gauss-Seidel 法････････････････････････・・・・・・・・・・・・・・･･･ 1 河野敏行（岡山理科大学），仁木滉（岡山理科大学），小武守恒（岡山理科大学），森本宗典（(株)日本情報管理システム） 15:40-16:05 非 H-行列に対する反復法の適用･･・・・・・・････・・･････・・・・・・・・・・・・・・・・・・・ 7 戸村健作(岡山理科大学)，仁木滉（岡山理科大学），森本宗典（(株)日本情報管理システム），榊原道夫（岡山理科大学） 16:05-16:30 3 次元力学空間における形態形成の最適化問題･････････・・・・・・・・・・・・・・・・ 13 有尾一郎(広島大学大学院工学研究科社会環境システム専攻)，石井崇之(日本総合研究所 (研究当時:広島大学大学院工学研究科構造工学専攻)) 16:30-16:55 Bilinear スキームの特性に関する考察と３次元移流方程式への応用・・・・・・・・ 19 松田洋介（熊本大学大学院自然科学研究科），畑上到（熊本大学大学院自然科学研究科）， 2 日目午前（5 月 7 日水曜日）第 2 セッション（9:00-10:15） 9:00-9:25 代用電荷法による標準スリット領域から多重連結領域への数値等角写像の方法・25 岡野大(愛媛大学工学部)，寺薗正俊(愛媛大学工学部)，緒方秀教(愛媛大学工学部)，天野要(愛媛大学工学部) 9:25-9:50 領域分割による SOR 法の収束性に関する考察･･･････････････・・・・・・・・・・・・・・ 31 畑上到（熊本大学大学院自然科学研究科），坂谷金哉（熊本大学大学院自然科学研究科）.

(2) 9:50-10:15 １次元カオス写像を用いたモンテカルロ積分･･･････････････・・・・・・・・・・・・・・・ 37 矢野忠（愛媛大学工学部），和田武（愛媛大学総合情報処理センター），矢野浩一（統計数理研究所），江沢康生（愛媛大学理学部） 10:15-10:30. ------ 休憩(15 分) ------. 第 3 セッション（10:30-12:10） 10:30-10:55 CG 画像の特徴を考慮した可逆画像圧縮法について･････････・・・・・・・・・・・・・ 41 井上啓 (山口東京理科大学基礎工学部)，黒塚豊(アルファシステムズ)，古市茂(山口東京理科大学基礎工学部) 10:55-11:20 二重量子ドットにおけるデコヒーレンスと散逸･･･････････････・・・・・・・・・・・・・・ 47 中村賢（山口東京理科大学基礎工学部） 11:20-11:45 高頻度金融時系列の特徴について････････････････････････・・・・・・・・・・・・・・ 53 田中美栄子（宮崎大学工学部） 11:45-12:10 A Parallel Algorithm for the Linear Complementarily Problem with an M-matrix・ 59 李磊(法政大学工学部) -------------------------------------------------------------------------------昼食 -------------------------------------------------------------------------------13:00 観光集合場所：稲盛会館前，ホテルウェルビュー前行先：桜島溶岩道路経由，桜島火山口，島津邸跡（磯庭園）予定 18:00 懇親会会場：稲盛会館.

(3) 3 日目午前（5 月 8 日木曜日）第 4 セッション（9:00-10:40） 9:00-9:25 A Markov modulated random walk with exponential upward ･･･・・・・・・・・・・・・・・・ 65 and exponential downward jumps 高田寛之(山口東京理科大学基礎工学部電子情報工学科） 9:25-9:50 ある種の任意次数微分方程式の解の微小摂動・緩和･･････・・・・・・・・・・・・・・・・ 71 都田艶子（大阪大学情報科学研究科） 9:50-10:15 座標仮定による有限要素法を用いた構造物の幾何学的非線形解析・・・・・・・・ 75 本間俊雄(鹿児島大学工学部)，大岡晃子(鹿児島大学理工学研究科)，安宅信行(昭和女子大学大学院生活機構学専攻) 10:15-10:40 Nekrasov 方程式の正値解の大域的な一意性に対する数値的検証法・・・・・・・ 81 小林健太（九州大学数理学研究院） 10:40-10:55. ------ 休憩(15 分) ------. 第 5 セッション（10:55-13:00） 10:55-11:20 移流拡散方程式の近似方程式に対する厳密解を用いた数値解法の提案・・・・ 85 榊原道夫（岡山理科大学），岡本直孝（岡山理科大学），仁木滉（岡山理科大学） 11:20-11:45 準残差の最小化によって加速パラメータを決める積型 BiCG 法・・・・・・・・・・・・・ 89 藤野清次 (九州大学情報基盤センター) 11:45-12:10 安定化近似逆行列法の改良･･･････････････････････････・・・・・・・・・・・・・・・・・ 95 池田優介 (九州大学大学院システム情報科学府) ，藤野清次 (九州大学情報基盤センター) 12:10-12:35 Family of cubic spiral transition curves ････････・・・・・・・・・・・・・・・・・・・・・・・・・ 101 Zulfiqar Habib and Manabu Sakai（Kagoshima University） 12:35-13:00 3 次元周期的領域における Stokes 流問題に対する基本解法･･･・・・・・・・・・・・ 103 緒方秀教（愛媛大学工学部情報工学科），安田雄二（愛媛大学大学院理工学研究科） 13:00. 閉会の挨拶.

(4) (I + R) 型前処理付 Gauss-Seidel 法河野敏行，仁木滉，小武守恒（岡山理科大学）森本宗典（ (株) 日本情報管理システム）. 1. はじめに線型方程式 Ax = b に前処理を施した方程式. P Ax = P b,. (1.1). について考える. ここで , A ∈ Rn×n は既約優対角 L-行列, P は正則な前処理行列, x, b はベクトルである. 本論文では , 一般性を失わないことから議論を簡単にするために A を A = I − L − U と置く. ただし , I は単位行列, −L と −U はそれぞれ A の狭義下三角, 狭義上三角行列である. 1991 年, Gunawardena らは , P = I + S を用いて , Gauss-Seidel 法を適用する修正反復法を提案した [2]. ただし , S は . S = (sij ) =. −aii+1 , 1 ≤ i ≤ n − 1, 0, otherwise,. である. さらに , 河野らは (I + S) 型の改良としてパラメータ α を用いた (I + αS) 型を提案し , α の推定値を導出している [3]. また, 小武守らは (I + S) 型の特徴を生かした (I + Smax ) 型の前処理を提案した [4]. ただし , Smax は . Smax =. (sm ij ). =. −aiki , 1 ≤ i ≤ n − 1, 0, otherwise,. ki = min{j | |aij | = maxi+1≤k≤n |aik |} である. これとは別に , A の狭義上三角行列 U を用いた適応的反復法 (I + U) 型が薄井らによって提案されている [8]. そして , 小武守らは (I + U) 型のかわりに (I + βU) 型を用いる前処理を提案している [5, 7] ただし , β は正の実数である. 上述の前処理法は , 与えられた係数行列の条件数を改善すると同時に , Gauss-Seidel 法の反復行列のスペクトル半径を減少させる効果がある. しかしながら , 行列 A の第 n 行に対して前処理の効果がない. 本発表では , 係数行列 A の第 n 行に前処理効果を与えるための前処理行列を提案し , 収束定理を証明する. 数値例にて，この前処理行列と上述の前処理行列とを併用し，収束率のより良い改善を試みる.. 1.

(5) 2. 新しい前処理行列の提案. n 行目に対する前処理効果として前処理行列 P = I + R, を提案する. ただし , R は . R = (rnj ) =. −anj , 1 ≤ j ≤ n − 1, 0, otherwise,. である. このとき, 前処理化係数行列 (I + R)A = AR の要素 aR ij は. aR ij. (2.2).    . =  . aij , anj −. n−1 k=1. 1 ≤ i ≤ n − 1, 1 ≤ j ≤ n, ank akj ,. 1 ≤ j ≤ n,. となる. ここで , AR を AR = MR − NR と分離する. ただし , MR は AR の下三角行列, NR は狭義上三角行列である. この分離を Gauss-Seidel 分離とよぶ . MR , NR はそれぞれ. MR = I − L + R − RL − RU, NR = MR − AR = U である.. n−1 k=1. ank akn = 1 のとき MR−1 が存在し , AR に対する Gauss-Seidel 反復行列 MR−1 NR. が定義できる. この論文では , n−1 . (2.3). k=1. ank akn = 1,. を仮定する. 次節で A = M − N, AR = MR − NR に対する比較定理を証明する. ただし , A = M − N, AR = MR − NR は Gauss-Seidel 分離である.. 3. 比較定理. まず , はじめに Gauss-Seidel 分離 A = M − N, AR = MR − NR が共に収束分離であることを示す. 補題 3.1 A = I −L−U を既約優対角 L-行列と置く. そのとき, A = M −N, AR = MR −NR は共に収束分離である.. 2.

(6) 証明まず , A = M − N が収束分離であることを示す. 仮定から , A が M-行列, すなわち A−1 > O を得る (Varga [9; Corollary 3.20]). 従って , A = M − N は収束分離である (Varga [9; Theorem 3.19]). 次に , AR = MR − NR が収束分離であることを示す. 仮定と (I + R) ≥ O より,. (I + R)Ax = AR x ≥ Ax ≥ 0 を満たす x ≥ 0 が存在する. したがって AR は単調行列である. すなわち A−1 R ≥ O を満たす (Axelsson [1; Lemma 6.1]). D, E をそれぞれ RU の対角行列, 狭義下三角行列と置いたとき,. MR−1 = [I + (I − D)−1 (L − R + RL + E) +{(I − D)−1 (L − R + RL + E)}2 + · · · +{(I − D)−1 (L − R + RL + E)}n−1 ](I − D)−1. (3.1). を得る. L ≥ R ≥ O より, L − R + RL + E ≥ O が成立する. さらに 0 ≤. n−1 . ank akn < 1 か. k=1. ら , (I − D)−1 ≥ O が成立し , MR−1 ≥ O を得る. したがって , AR = MR − NR は正則分離であり, 収束分離である (Varga [9; Theorem 3.29]). 補題 3.2 A を既約優対角 L-行列と置く. さらに A = M − N を Gauss-Seidel 正則分離と置く. そのとき, T z = ρ(T )z を満たす正のベクトル z が存在する. ここで , T = M −1 N である. 証明 T は以下のように表せる. . T =. 0 T0 0 T1. . T1 は (n−1)×(n−1) 行列である. 仮定から , T1 は既約で非負行列であるから Perron-Frobenius の定理より, ρ(T1 ) = 0, T1 y = ρ(T1 )yを満たす y > 0 が存在する. z を . z=. 1 Ty ρ(T1 ) 0. y.  . と置く. . Tz =. T0 y T1 y. . . , ρ(T1 )z =. T0 y ρ(T1 )y. . . =. T0 y T1 y. . ,. より, T z = ρ(T1 )z が成立し , ρ(T1 ) = ρ(T ) となり, T z = ρ(T )z を満たす正のベクトル z が証明された.. 3.

(7) 補題 3.3 [2] T を非負行列とする. そのとき,. (a) ある非負ベクトル z (z = 0) に対して αz ≤ T z を満たすとき, α ≤ ρ(T ). (b) ある正のベクトル z に対して T z ≤ βz を満たすとき, ρ(T ) ≤ β, さらに T が既約で , ある非負ベクトル z に対して ,. 0 = αz ≤ T z ≤ βz のとき, α ≤ ρ(T ) ≤ β であり, z は正のベクトルである. 定理 3.4 A を既約優対角 L-行列と置く. さらに A = M − N, (I + R)A = MR − NR を共に Gauss-Seidel 正則分離と置く. そのとき,. ρ(MR−1 NR ) ≤ ρ(M −1 N) < 1 が成立する. 証明仮定と補題 3.2 より次式が成立する.. 1. Mz =. (3.2). ρ(M −1 N). Nz ≥ 0. A = M − N は収束分離であるから次式を得る. Az = (M − N)z = M(I − M −1 N)z =. 1 − ρ(M −1 N) Nz ≥ 0. ρ(M −1 N). したがって, R ≥ O と N = NR = U より. MR (M −1 − MR−1 )Mz = {(MR − N) − (M − N)}z = {(I + R)A − A}z (3.3). = RAz ≥ 0. が成立する. したがって, (3.2) から次の不等式が成立する.. (3.4). MR (M −1 − MR−1 )Nz ≥ 0.. MR−1 ≥ O より (M −1 − MR−1 )Nz ≥ 0, となり,. (M −1 − MR−1 )Nz = M −1 Nz − MR−1 NR z = ρ(M −1 N)z − MR NR z ≥ 0 が成立する. 補題 3.3 (b) と ρ(M −1 N) < 1 より, ρ(MR−1 NR ) ≤ ρ(M −1 N) < 1 が成立する. 発表にて，数値例を示す．. 4.

(8) 参考文献 [1] O. Axelsson, Iterative Solution Methods, CAMBRIDGE UNIVERSITY PRESS, 1994. [2] A. D. Gunawardena, S. K.Jain and L. Snyder, Modified Iterative Methods for Consistent Linear Systems, Linear Algebra Appl. 154-156(1991)123-143. [3] T. Kohno, H. Kotakemori, H. Niki, Improving the Modified Gauss-Seidel Method for Z-matrices, Linear Algebra Appl., 267(1997)113-123. [4] H. Kotakemori, K. Harada, M. Morimoto and H. Niki, A comparison theorem for the iterative method with the preconditioner (I + Smax ), J. Comput. Appl. Math., 145(2002)373-378. [5] H. Kotakemori, H. Niki and N. Okamoto, Accelerated iterative method for Z-matrices, J. Comput. Appl. Math., 75(1996)87-97. [6] M. Morimoto, K. Harada, M. Sakakihara, H. Sawami, The Gauss-Seidel iterative method with the preconditioner (I + S + Sm ), Japan Journal of Industrial and Applied Mathematics,(To appear). [7] 薄井正孝, 仁木滉, 小武守恒, 河野敏行, (I + βU) 型前処理付 Gauss-Seidel 法の収束定理, 日本応用数理学会論文誌, Vol.6, No.4,(1996) 307-316. [8] M. Usui, H. Niki, N. Okamoto, Adaptive Gauss Seidel method for linear systems, Intern. J. Computer Math., 51(1994) 119-125. [9] R.S. Varga, Matrix Iterative Analysis, Springer, 2000.. 5.

(9) 6.

(10) 非Ｈ行列に対する反復法の適用戸村健作 ∗ ，森本宗典 ∗. 1. ，榊原道夫. ∗∗∗. ，仁木滉. ∗∗∗. 岡山理科大学大学院総合情報研究科 ∗∗. ∗∗∗. ∗∗. （株）日本情報管理システム. 岡山理科大学総合情報学部情報科学科. はじめに線型方程式. Ax = b. (1). を反復法で解くことを考える．ただし，A ∈ Rn n は正則な行列，x，b はそれぞれ未知，既知のベクトルである．A の分離を A = M − N とする．ただし，M は正則行列である．このとき，(1) の反復公式は. x(k+1) = M −1 N x(k) + M −1 b, (k = 1, 2, · · ·) ,. (2). となる．係数行列 A が H 行列のとき，A の反復行列 M −1 N のスペクトル半径が ρ(M −1 N ) < 1 となることが知られている [1][5]．Gauss-Seidel 法では M ，N はそれぞれ下三角行列，狭義上三角行列である．Gauss-Seidel 法による反復公式を考えた場合，A が非 H 行列のとき，収束性が保証されないために，反復法は適用できない．そのために Gauss の消去法が用いられているが，計算量と大型行列のときに発生する丸め誤差の問題がある．そこで我々は，係数行列が非 H 行列のとき，前処理法を用いることで H 行列に変換する方法を提案する．従来の前処理法には，Gunawardena らによる P = I + S 型前処理 [2]，小武守らによる P = I + Smax 型前処理等 [4] があり，これらは H 行列に対する加速処理として用いられてきた．今回，我々は係数行列 A が非 H 行列のとき，(1) 式に左から正則な前処理行列 P で前処理化で，得られた方程式 P Ax = P b の前処理化係数行列 P A を H 行列化する前処理法を提案する．. 2. 定義と定理. 定義 2.1 A = (aij ) ∈ Rn×n に対して，aij ≤ 0, i 6= j を満たす行列を Z 行列という．また， aii > 0 (1 ≥ i ≥ n) を満たす Z 行列を L 行列という．また，A が正則で A−1 ≥ 0 のとき，この行列を M 行列という．. 7.

(11) 定義 2.2 A = (aij ) ∈ Rn×n に対して, (. αij =. |aij |, for i = j −|aij |, for i 6= j. とおく行列を比較行列と呼び，< A >= (αij ) で表す． < A > が M 行列のとき，つまり < A >−1 ≥ O のとき A を H 行列という．定義 2.3 [3] 行列 A を A = M − N と分離する．ただし，M は正則である．そのとき，. (a) < M > −|N | が M 行列のとき，H 分離， (b) < A >=< M > −|N | のとき，H 互換分離, と定義する. 定理 2.4 [3] 分離 A = M − N が，H 分離のとき, A と M は H 行列であり， ρ(M −1 N ) ≤ ρ(< M >−1 |N |) < 1 が成立する. 定義 2.5. A = (aij ) ∈ Rn×n に対して |aii | ≥. n X. |aij |. (3). j=1 j6=i. 成立するとき, 優対角行列という. (3) 式ですべての i について, 狭義の不等式が成立するとき, 狭義優対角行列という. (3) 式において少なくとも１つの i に対して狭義の不等式が成立するとき, 既約優対角行列という. 定義 2.6 行列 A の第 i 行の優対角度 ti を， n X. ti =. |aij |. j=1 j6=i. |aii |. (i = 1, 2, · · · , n). と定義する．. 定理 2.7 n 次正方行列 A が狭義優対角行列，または既約優対角行列のとき，A は H 行列である．. 8.

(12) 3. 非 H 行列の H 行列化非 H 行列 A に P で前処理を施した行列 P A の H 行列化の方法を提案する．. 前処理法１前処理行列を P = P (k) · · · P (3) P (2) P (1) とし，A(k) = P A とおく．ただし，A(0) = A である． P (k) は，   1, j=i     (k−1)  aimi (k) (k−1) P (k) = (pij ) = − (k−1) , j= 6 i, ti >1   ami mi     0, otherwise (k−1). mi = min{j : |aij. | = 1≤l≤n max |ail |} l6=i. (k−1). を満たす行列である．ただし，ti は k − 1 回前処理を施した係数行列 A(k−1) の優対角度である．P (k) は A(k−1) の各行の絶対値最大の要素を 0 にする前処理効果がある．P (k) の前処理を繰り返すことで非優対角行列を優対角行列に変換する．. · · · P (3) P (2) P (1) Ax = · · · P (3) P (2) P (1) b. 前処理法２ (k). (3). (k−1). (2). (1). · · · PU PL PU とし，A(k) = P A とおく．ただし，A(0) = A 前処理行列を P = PL PU (k) (k) である．ここで，PL ，PU は，それぞれ以下に示す行列である．   1     (k−1)  a i (k)L (k) PL = (pij ) =  − im (k−1)  ami mi     0 (k−1). mi = min{j : |aij. (k−1). (k). (k−1). j < i, ti. > 1,. otherwise. (k−1). | = max |ail 1≤l≤i−1.   1     (k−1)  a i (k) (k)U PU = (pij ) = − im (k−1)   ami mi     0. mi = min{j : |aij. j = i,. |}. j = i, (k−1). j > i, ti. > 1,. otherwise. (k−1). | = max |ail i+1≤l≤n. |}. (k). PL ，PU はそれぞれ A(k−1) の狭義下三角，狭義上三角行列の各行の絶対値最大の要素を 0 に (k) (k) する．PL ，PU の前処理を交互に繰り返すことで非優対角行列を優対角行列に変換する． (3). (2). (1). (3). (2). (1). · · · PU PL PU Ax = · · · PU PL PU b. 9.

(13) 前処理法１，２により得られた線型方程式. P Ax = P b. (4). の係数行列は優対角行列なので定理 2.7 より H 行列である．(1) と (4) は等価であるから，(4) を Gauss-Seidel 法で解くことを考える．P A = Mp − Np を P A と分離する．このとき，Mp を下三角行列，Np を狭義上三角行列とする．P A は H 行列であるから，< P A >=< Mp > − < Np > は H 分離であり，H 互換分離であることがわかる．よって，定理 2.4 より. ρ((Mp )−1 Np ) ≤ ρ(< Mp >−1 |Np |) < 1 となり A(k) に対する Gauss-Seidel 法は収束する．. 4. 数値例. 我々が提案した方法の妥当性を数値例で確かめる． 2 次元ラプラス方程式. ∇2 φ = 0 in Ω = (0, 1) × (0, 1) を境界要素法による離散化方程式 Aφ = b をテストする. 係数行列 A は非 H 行列なので, H 行列に変換し, Gauss-Seidel 法で得られた結果を Gauss の消去法と比較する. 収束判定条件は ||φ(k+1) − φ(k) ||2 /||φ(k+1) ||2 ≤ 10−12 とする. 境界要素解析は一定要素を用い，境界条件を. φ(0, y) = 1.0, φ(x, 0) = φ(1, y) = φ(x, 1) = 0.0 とする. このとき, 400 から 2000 までの要素に対する結果をを表 1 に示す. 表 1.H 行列化反復法の計算結果要素数. 平均優対角度. 400 800 1200 1600 2000. 37.19 67.81 96.56 124.21 151.10. Gauss 消去法実行時間 [s] 0.192 1.725 5.727 13.444 25.964. 要素数. 前処理回数. 400 800 1200 1600 2000. 299 783 1255 1815 2394. 前処理回数. 146 374 596 908 1258. 前処理法１ + Gauss-Seidel 法前処理時間 [s] 反復反復時間 [s]. 0.177 1.389 4.701 12.708 26.561. 18 18 18 18 18. 前処理法２ + Gauss-Seidel 法前処理時間 [s] 反復反復時間 [s]. 0.109 1.082 3.906 9.475 18.762. 19 19 19 19 19. 0.012 0.048 0.109 0.193 0.301. 0.011 0.048 0.103 0.183 0.333. 合計時間 [s]. 0.188 1.438 4.804 12.891 26.895. 合計時間 [s]. 0.122 1.131 4.016 9.669 19.063. (使用 PC CPU Pentium4 2.4MHz, memory 512MB). 10.

(14) 5. 終わりに. 本研究では，非 H 行列を係数行列にもつ線型方程式を，２つの前処理法で H 行列に変換し反復法の適用可能化を試みた．その結果，前処理法１，２共に H 行列化に成功し，Gauss の消去法よりも速く近似解を得ることができた．前処理法１の特徴は，少ない前処理回数で H 行列に変換できることである．前処理法２は前処理法１よりも前処理回数は増えるが，CPU タイムが減少した．. 参考文献 [1] A.Berman and R.J.Plemmons, Nonnegative Matrices in the Mathmatical Sciences ,SIAM,philadelphia(1994). [2] A.D.Gunawardena, S.K.Jain, and L.Snyder, Modified Iterative Methodsfor Consistent Linear Systems, LinearAlgebraanditsApplications., 154-156(1991)123-143. [3] A.Frommer, D.B.Szyld, H-splittings and two-stage iterative methods,Numerische Mathematik.63 (1992) 345-356. [4] H.Kotakemori, K.Harada, M.Morimoto, and H.Niki, A comparison theorem for the iterative method with the preconditioiner(I +Smax ), J.Comput.Appl.M ath.,145(2002)373-378. [5] R.S.Varga, Matrix Iterative Analysis 2nd edition, Springer(2000).. 11.

(15) 12.

(16) Optimization Lay-out Problem in 3-dimensional Finite Element Space 1 I. ARIO∗1 and T. ISHII∗2 ∗1. Dept. of Civil and Environmental Engineering, Hiroshima University, 1-4-1 Kagamiyama, Higashi-Hiroshima, 739-8527, Japan ∗2 The Japan Research Institute, Limited, 16 Ichiban-cho, Chiyoda-ku, Tokyo. Abstract This paper presents an iterative structural design technique that uses local control conditions and the stiffness of individual members. A lattice truss is modeled as an assemblages of a number of repeating unit cells to undergo finite element analysis in order to feed-back the results for the stiffness of each member. Further iteration enable for the identification of a topology with the use of the local stress conditions without mathematical optimization. The analysis of widely used shape-layout problems has proposed technique. Key Words: nonlinear equilibrium equation, structural optimization, shape-layout problem, iterative stiffness control. 1. Introduction. The work of Bendsøe and Kikuchi (1988) on continuum body topology optimization using the homogeneous material method, and recent advances in computing technology have enabled for a number of significant researches in the field of structural design. For example [3, 4] proposed the repeating method and Honma and Tosaka[5] presented condition equations for structural control and performed an analysis of structural design using repeated calculations. Xie & Steven[6] developed Evolutionary Structural Optimization (ESO), which based on repeated Finite Element Method calculation and shape control. However problems with this technique include the selection of unnecessary members and the replacement of unnecessary members after the process has finished. Ohmori et. al [7] proposed deletion standards using stress contours and local areas, and suggested an expanded ESO method with one-sided deletion. In this paper, we propose a structural design method using simple local rules of nonmathematical optimization similar to the ESO method. It is a nonlinear repeating system with a feedback response system for design variables (stiffness). The identification of an optimum design is difficult with structural optimization methods such as the mathematical optimization and sensitivity analysis methods due to control conditions such as stress and weight. Even when such conditions are satisfied, the result is not always a nonlinear global optimum solution. The proposed system uses an expanded design area in which elements can be revived, a process that is not possible with the original ESO method. In addition, the stress response to design variables can be feed back into the system.. 2. Structural Design Method. In this section, the concept of the design method based on local stresses member is introduced. 1. 13. The part of this manuscript have been presented in the citation [1].

(17) 2.1. Concept of layout problem. We define a truss design area Ω containing finite design variables as . T. x = · · · , x(m) , · · ·. ∈ RM in Ω. (1). The nonlinear equation is defined as F (u, p, x) = 0. (2). Now u ∈ RN is the displacement vector, and p ∈ R is the loading parameter. Finally, we search for solutions (u, p, x) that satisfy the discrete equilibrium equation(2). If we locally consider linear problem near to the equilibrium point in equation (2),the incremental equation is ∂F ∂F + =0 Ju p + x (3) ∂p ∂x Where • denotes incremental variables. Let J(u, p, x) =. ∂F ∂u. be the tangent stiffness matrix. However, then this objection to this method is to structural optimization, for example minimum weight maximum stiffness [8] is described as M . x(m) → min. (4). m=1. s.t. xmin ≤ x(m) ≤ xmax , σmin ≤ σ (m) ≤ σmax , umin ≤ u(m) ≤ umax , ∀ m = 1, · · · , M.. (5). The main aim of this technique is not the identification of a mathematical optimum, but to identify the local stress control conditions of each member. We identify a topology calculating by the stiffness of member. If the aim of the structural optimization is to minimize weight, the main problem is to balance the entire structure system. In this research, stress share control design variable, consequently the proposed method identifies a rational structure which satisfies the stress condition.. 2.2. Structural design using to a section feedback control. When change of stiffness alternation is limited part of the structure, the structure can be analyzed with a higher degree of accuracy using the direct method. However the use of an iteration method is more applicable for the modification of an entire structure. Modification of the stiffness matrix is defined as feedback which provides each member with a stress response. We control the stiffness, to satisfy the equilibrium equation. Let us establish an initial design variable x(0) and an equilibrium point (u, p, x(0) ) for weight or displacement control. The stress in each member is defined as σmin ≤ σ (m) ≤ σmax , . . σ (ν) = W u(ν) ,. m = 1, · · · , M ν = 0, 1, · · ·. (6) (7). 14.

(18) . It is the displacement function of each node. Here ν is defined as the number of iterations, (m) (m) and σ (ν) = (· · · , σ(ν) , · · ·)T , is defined as u(ν) = (· · · , u(ν) , · · ·)T . Furthermore, the present design variable (stiffness) x(ν) is resolved into x(ν+1) . We then define the resolve rate γ. A new design variable is also resolved into . . x(ν+1) = F γ, σ (ν) , . . = F γ, W u(ν). . ,. ν = 0, 1, · · ·. (8). We ensure that the stiffness matrix is also renewed. This result enable us to identify solutions (u(ν+1) , p, x(ν+1) ) of the equilibrium equation (ν+1) = J(u, p, x(ν+1) )u. ∂F p ∂p. (9). 2. We substitute this solution into equation (7) and perform a series of repeated calculation, such that the convergence stress and displacement conditions are sat. In other words, the nodal displacement and design variable are expressed as . . u(ν) = F p, x(ν) , . . x(ν+1) = F γ, W(u(ν) ) ,. ν = 0, 1, · · ·. (10). by repeated calculation of the equilibrium equation. If γ, p, W is included, multi-dimensional and multi-feed back of the nonlinear repeat calculation is expressed as x(ν+1) = F (x(ν) ) = F (F (· · · F (x(0) ))), = F ν (x(0) ), ν = 0, 1, · · · ,. (11). . We adopt a limited stress condition and variable based on local rules.. 3. Examples. 3.1. Problem of the coat hook. The design of a simple coat hook subjected to vertical load p at a free point far from a wall, is a widely known structural design problem. In this case, we try as the bench model of the layout formation using the proposed analysis method. The lattice structure which is shown in Fig.1 is modeled using 16 × 16 elements, and we establish a fixed point whose aspect ratio is 1.5. This fixed point is third node to the inside from the top and bottom edges of the design area. All initial member stiffness which is design variables are (EA0 = 1). The maximum stiffness is assumed to be EAmax = 10, and the load parameter initiates stiffness control under p/EA0 = 0.0001. The stress with respect to time is shown in Fig.2. In this figures member thickness is a magnitude of stress. Horizontal members at the center and corners adopt a zero stress state in the initial step. A void is observed in Fig.2(c), after which reinforcement members can be seen. In particular, final truss structure was formed of the compression member between two tension members near to the bearing point. Finally, an ideal frame structure which is made from up members with the maximum set condition member stiffness can be seen in Fig.2. The contents in this section will be presented on the day. 2. 15. Now this solution is defined as displacement vector u(ν+1) for weight control..

(19) 10. 7.5. 5. 2.5. 0 0. 2.5. 5. 7.5. 10. Figure 1: Lattice truss of 16 × 16 cells. (a) step 0. (b) step 2. (c) step 4. (d) step 6. (e) step 8. (f) step 9. Figure 2: Process of optimum formation of 16 × 16 cells. 16.

(20) yy ;; ; y ;; ;; yy y ;; yy yy; ;; yy yy ;; Figure 3: Michell’s problem in 3-dimensional space. (a) forming at step 4. side view at step 4. (b) forming at step 10. side view at step 10. (c) forming at step 30. side view at step 30. Figure 4: Process of the arch structure forming. 17.

(21) 4. Conclusion. We have shown a process for the development of truss structures using the forces in members as a force flow in the design area, and fed back this information as a structural variable without mathematical optimization or the sensitivity analysis method. Especially, we could show possibility that we can make structural design or shape (structural beauty or smart structures) by computing analysis method in this research. For example efficacy of structural design and change for small or light. Further research will focus on dynamic solutions for optimum material layouts and slim structures.. References [1] Ario, I. and Ishii, T. : Structural Design using Repeated Stiffness Control, Proceedings of China-Japan-Korea Joint Symposium on Optimization of Structural and Mechanical Systems, 2, (2002), 821-826. [2] Bendsøe, M. P. and Kikuchi, N. : Generating Optimal Topologies in Structural Design using a Homogenization Method, Computer Methods in Applied Mechanics and Engineering, 71, (1988), 197-224. [3] Suzuki, K. and Kikuchi, N. : A homogenization method for shape and topology optimization, Computer Methods in Applied Mechanics and Engineering, 93, (1991), 291-318. [4] Schwarz, S., Maute, K. and Ramm, E. : Topology and shape optimization for elastoplastic structural response, Comput. Methods Appl. Mech. Engrg. 190, (2001), 21352155. [5] Honma, T. and Tosaka, N. : Ahape Analysis and Optimization of Structure using an Iterative Caluculation Method, J. Struct. Constr. Eng., AIJ (in Japanese), 502, (1997), 69-76. [6] Xie, Y.M. and Steven, G.P. ：Evolutionary Structural Optimization, Springer, 1997. [7] Ohmori, H., Cui, C., and Suzuki, N.：Creation of Structural Form by Extended ESO Method Trough Usage of Contour Lines, Structural Engineering Symposium, AIJ (in Japanese), 47B, (2001), 7-14. [8] M. P. Bendsøe : Optimization of structural topology, shape and material, Berlin, Heidelberg, New York : Springer, 1995.. 18.

(22) Bilinear スキームの特性に関する考察と３次元移流方程式への応用松田洋介1 ，畑上到（熊本大学大学院自然科学研究科）. 1. 緒言双曲型の偏微分方程式の代表的なものである移流方程式を差分化により, 数値的に解こうとすると, スキームの特性が顕著に反映され, 様々な影響が出てくる. 例えば空間１次精度のスキームでの物理量の散逸や, 空間２次精度のスキームで現れる分散などによる初期状態からの形状の変化などがその代表的な例である. このような現象が生じる原因として，一般的に用いられているスキームが, 数値流束の評価にあたり各座標軸の方向毎に独立にしか考慮しておらず, 多次元的な項（例えば斜め方向についての項）を評価していないためであると考えられる. このことは, スキームの誤差という点でも, 交差方向成分に対して大きな数値誤差が入るという問題も引き起こす. したがって，このような問題に対処していく上では, あらゆる方向に対して同時に考慮に入れたスキームが必要になってくる. 本研究では, 時間発展に対して初期条件の形状を維持することが期待されるスキームとして Shubin らにより提案された Bilinear スキーム 1) に着目し, その特性について詳細に解析した．Bilinear スキームとは, 数値流束を計算する際に, 各座標軸毎に分割して計算を行う手法をとるのではなく, 各座標軸を互いに結合させた形で, 多次元的に評価を行い計算するスキームであることから, そのような名称で呼ばれる. さらに，本研究では, この Bilinear スキームの特性を十分活かして３次元系への拡張を試みた.. 2. 解析方法 Bilinear スキームについて２次元の場合について簡単に説明する．線型の２次元移流方程式の初期値問題   ∂u ∂(c1 u) ∂(c2 u) + + = 0 (t > 0, −∞ < x, y < ∞, c1 , c2 : 定数) (1) ∂x ∂y  ∂t u(x, y, 0) = ϕ(x, y) (−∞ < x, y < ∞) を，以下のように離散化する． n un+1 i,j = ui,j. ∆t − ∆x (c1 ui+1/2,j − c1 ui−1/2,j ) ∆t − ∆y (c2 ui,j+1/2 − c2 ui,j−1/2 ). .. (2). Bilinear スキームでは, 格子点間の数値流束，ui±1/2,j , ui,j±1/2 を, リーマン解法 (Riemann Solver) の考え方により，. ui±1/2,j = ui,j±1/2 =. uL i±1/2,j. if. c1i±1/2,j ≥ 0. uR i±1/2,j. if. c1i±1/2,j < 0. uL i,j±1/2. if. c2i,j±1/2 ≥ 0. if. c2i,j±1/2 < 0. uR i,j±1/2. (3) (4). として，これらを多次元的に決めることにより, 適合性を満たし, かつ多次元的に空間２次精度となる Bilinear スキームを構築する. まずここでは，x 方向の微分について，c1 が正のとき, つまり uL i+1/2,j について考え. L る. 以下に uL i+1/2,j を評価する考え方を示す. 図 1 のような領域 ABCDEF で ui+1/2,j を考える．この図. は, 格子点 i, j を中心とする xy 平面上の領域 Bi,j において, 時間 tn+1 の i + 1/2 の領域（図では直線 CF. 上）に, 時間 tn の, どの領域の物理量 u が伝わるかを, 空間-時間的に示している. また時間 tn における領 1. 現在の所属は日産自動車（株）. 1. 19.

(23) 域 ABED に存在する u の情報が, 時間 tn+1 の i + 1/2 の領域に伝わることを示しており, 時間 tn から微少時間 ∆t 後の tn+1 までには ABCDEF 内の u の情報が影響を及ぼすことを意味している.. F. −c1 ∆t, ∆y , tn ) ( ∆x 2 2 Bi,j. y. ( ∆x , ∆y , tn+1 ) 2 2. ( ∆x , ∆y , tn ) 2 2 D. E C. , − ∆y , tn+1 ) ( ∆x 2 2. B. ( ∆x , − ∆y , tn ) 2 2. t. A ∆y n ∆x ( 2 −c1 ∆t, − 2 , t ). x. 図 1: x 方向についての概念図ここで，(1) を図 1 の領域 ABCDEF で積分することにより，次の式を得る. ut + c1 ux + (c2 u)y dxdydt = 0. (5). ABCDEF. ここで，c2 u をまとめているのは，c2 が定数でない場合にも適用するようにするためである．これを項毎に積分すると,. c1 ∆t∆yuL i+1/2,j. = c1 u dydt BCEF = u dxdy + ABDE. ABC. (6) (c2 u) dxdt −. この式において, 右辺第一項は中間点の値 uM を使って, u dxdy = c1 ∆t∆yuM ABDE. (∆x − c1 ∆t) ux,i,j 2 と表すことができる. これより式 (6) は次のように表すことができる. uM = ui,j +. uL i+1/2 = uM +. ∆t − (Γ − Γ+ ) 2∆y. DEF. (c2 u) dxdt. (7). (8). (9). 2. 20.

(24) ここで Γ− ,Γ+ はそれぞれ c2 u の△ ABC, △ DEF における項であり, 以下のようにして求める（ここでは. DEF 面を考える）. 今，c1 , c2 がともに正であるとするとき, 図 2 のような四面体を考えることで Γ+ の評価を行う. これは時間ステップ tn から tn+1 の間に三角形 DEF がどの領域から影響を受けるかを表した概念図であり, 三角形 DEF は四面体 DEFG の内部に存在する全ての u から影響を受ける. つまり四面体 DEFG の内部が影響を与える領域となるので, そこに存在する u で評価を行う. また, この Γ+ や Γ− を考えるこ ∂(c1 u) の差分式の中に，y 方向からくる u の情報を取り入れることになり，Bilinear スとは,(1) における ∂x キームが, 多方向から同時に情報をとり入れるという特徴が現れるところである.. ( ∆x , ∆y , tn+1 ) 2 2. F. −c1 ∆t, ∆y , tn ) ( ∆x 2 2. D. E. ( ∆x , ∆y , tn ) 2 2. Bi,j. y. t. G ∆y ∆x ( 2 −c1 ∆t, 2 −c2 ∆t, tn ). x. 図 2: △ DEF についての概念図. m(DEF ), m(DEG) をそれぞれ ∆DEF ,∆DEG の面積を表すとすると，中間点の値を用いることにより, 1 L uDEF = u dxdy = um m(DEG) DEG ∆y 2 1 とできる. ただしここでの um は (x, y) = ( ∆x 2 − 3 c1 ∆t, 2 − 3 c2 ∆t)（重心）での u の値で, 具体的には. um = ui,j. 1 ∆t + 3∆x−4c uxi,j 6 2 ∆t + 3∆y−2c uyi,j 6. (10). R と表される. 以上により求まる uL DEF ,uDEF で Γ を定義することができる.. Γ+ = c2 uDEF uDEF =. uL DEF uR DEF. if if. (11) c2 ≥ 0 c2 < 0. (12). となる．三角形 ABC についても j − 1 の領域について評価することにより同様に求めることができ, 同じ L R く uR i+1/2,j も i + 1 の領域での評価で相似的に求めることができる. また ui,j±1/2 , ui,j±1/2 についても同様. に求められる. 以上の議論は，c1 , c2 が定数でない場合でも, 線形化によって同様に進められるが, そのため. 3. 21.

(25) に，c1 , c2 についての誤差が加わることになる. また，３次元への拡張は，かなり複雑になるので省略するが，同様に構築できる．. 3. 計算モデルおよび条件本研究では, 一辺の長さが 6 の正方形領域について考え, 格子数が 501 × 501 である直交等間隔格子に分割して計算を行った. x, y 各方向の特性速度を c1 = c2 = 1.0 とし, 初期条件には, 次の式で表されるものを選択した.. u(x, y, 0) = exp[−60{(x − 3)2 + (y − 3)2 }]. (13). また境界条件については，c1 = c2 = 1.0 であることより, 進行方向に対して平行に, 斜めの値で境界条件を施すべきであるが, 今回の計算においては, 時間発展において, 境界付近ではなく, 初期の段階での形状の変化を確認するため, 各方向ごとにノイマン条件を与えた. 以上の条件の下で，Bilinear スキームの特性を考察するにあたり, 比較するスキームとして Lax-Wendroff のスキーム（以下略して，LW スキームと記す）を採用し，比較を行った. また３次元の場合にも，計算モデルとして，一辺が 1 の立方体領域を考え, 計算格子数が 101 × 101 × 101 の直交等間隔格子に分割し，x, y, z 各方向の特性速度を c1 = c2 = c3 = 1.0 とした．初期条件は以下の式で表されるものを選択した.. u(x, y, z, 0) = exp[−60{(x − 0.5)2 + (y − 0.5)2 + (z − 0.5)2 }]. (14). また，境界条件は, 全境界面でノイマン条件とした.. 4. 計算結果と考察 4.1 ２次元の移流方程式への Bilinear スキームの適用以上の計算条件での計算結果を以下に示す. まず, 図 3 に,Bilinear スキームと L-W スキームの両方について, 時間 T = 0.36 での u の形状の変化の様子を等高線で示した.. (a) Bilinear スキーム. (b) L-W スキーム. 図 3: 2 次元移流方程式における数値計算結果（等高線）. この結果から, L-W スキームが, 初期条件の形状から, 特に斜め方向に大きく変形しているのに対し,. Bilinear スキームは, 多少の変形は見られるものの, L-W スキームほど大きな変形は見られないことが分かった. また, スキームの特性を示す指標の一つとして,(15) で表される, 厳密解との誤差を各時間ステップ. 4. 22.

(26) で計算したものを図 4 に示す.. . E=. |ui,j − u(x, y, t)|. (15). i,j. L-W スキームの方が時間発展とともに，誤差が大きく増大していることがわかる．. Residual E of Bilinear and LaxWendroff 120 ’Bilinear_2D.dat’ ’LaxWendroff_2D.dat’ 100. Residual E. 80. 60. 40. 20. 0 0. 0.5. 1. 1.5. 2. 2.5. Time. 図 4: Bilinear スキームと L-W スキームの厳密解との誤差の時間的変化の様子. 4.2 Bilinear スキームと L-W スキームの誤差項の比較前節で示した計算結果で確認された Bilinear スキームの特性について, 両スキームをテーラー展開して，打ち切り誤差の持つ意味についての比較を行った. それぞれの主要な誤差項により，Bilinear スキームと L-W スキームでは, ut + (c1 u)x + (c2 u)y = 0 を解いているのではなく, 厳密には次のような式を解いていることになる.. ∂u ∂(c1 u) ∂(c2 u) + + = −αb ∂t ∂x ∂y ∂u ∂(c1 u) ∂(c2 u) + + = −αl ∂t ∂x ∂y. : Bilinear スキーム. (16). : L − W スキーム. (17). ここで,αb , αl はそれぞれ， 1 −{ 12 (∆t)c1 (∆x)2 − 14 (∆t)2 c21 (∆x) + 16 (∆t)3 c31 }uxxx (x, y, t). αb =. 1 −{ 12 (∆t)c2 (∆y)2 − 14 (∆t)2 c22 (∆y) + 16 (∆t)3 c32 }uyyy (x, y, t). αl. =. c1 c2 (∆t)2 uxy 1 1 −{ c31 (∆t)3 − c1 (∆t)(∆x)2 }uxxx 6 6 1 1 −{ c32 (∆t)3 − c2 (∆t)(∆y)2 }uyyy 6 6 1 1 − c21 c2 (∆t)3 uxxy (x, y, t) − c1 c22 (∆t)3 uxyy (x, y, t) 2 2. (18). (19). である．これから，L-W スキームには２階の微分項 uxy の誤差が主要な項として存在しており, この項が形状の変形に大きく影響していることが予想される．このため，この項を相殺するように項を加えて計算を行うと，変形は低く抑えられることがわかった．したがって, 元来そういった誤差項を持たない Bilinear スキームは, 誤差項を持つスキームより, 変形をすることなく初期条件の形状を維持することに長けていると考えられる.. 4.3 ３次元の移流方程式への Bilinear スキームの適用次に，Bilinear スキームを, ３次元の場合に適用した結果について示す. 図 5 には Bilinear スキームと,L-W スキームで計算した場合の,u の分布を, その時間における球の中心を通る yz 平面で切った面の等高線で示. 5. 23.

(27) し, 時間発展とともに並べた図で, 手前から, 時間 T = 0, 0.1, 0.2, 0.3 の状態を表している. この図より, 若干ではあるが，Bilinear スキームの方が形状を維持していることが分かる. さらに図 6 の Bilinear スキームで計算した場合と L-W で計算した場合の, 厳密解との誤差が時間発展に対して変化する様子を見ると,. Bilinear スキームの方が誤差を抑えていることが確認される.. (a) Bilinear スキーム. (b) L-W スキーム. 図 5: yz 平面で見る時間発展. Residual E of Bilinear and LaxWendroff 500 ’Bilinear_3D.dat’ ’LaxWendroff_3D.dat’ 450 400. Residual E. 350 300 250 200 150 100 50 0 0. 0.05. 0.1. 0.15 Time. 0.2. 0.25. 0.3. 図 6: 厳密解との誤差の時間的変化の様子. 5. 結論本研究では, ２次元あるいは３次元の移流方程式を計算する際に, 数値流束 ui±1/2,j , ui±1/2,j,k 等を２次元的, あるいは３次元的に評価するために, Shubin らによる Bilinear スキームに着目し, その特性について, １次元的な考え方の延長である, L-W スキームと比較することにより詳細に調べた. この Bilinear スキームでは, ２次元, あるいは３次元的に, 各座標軸成分を結合させて数値流束を評価することにより, 振動は若干生じるものの Lax-Wendroff のスキームより, 初期条件の形状を, より保持することができることがわかった. これは,uxy などの誤差項の影響によるものであり, この交差項を抑えることが重要であることがわかった.. 参考文献 1) J.B.Bell,C. N.Dawson and G. R.Shubin, J.Comput. Phys. 74,pp.1-24(1988). 6. 24.

(28)

(29) "!#$%&'( )+* , -/. 021 354 687 9 * : ; < ,2=2,>=2?A@2B5=DC>E2F ;>< ,2=DBG=DH>I2J2B5=DF (2 1. 1. 1. 1. ). K L>MON PRQRSUTWV yoxs RRXRRY[uRZ]vU\_W^[cZaR`RWbdgc dcgedcgqf+lmhj i+km lonq p R T cszRrR|RtRXRuR}vxw+Ry{zR |R XRym}Rx~s g c ¡U¢Rcs£¤cjfhiV k¥c§¦g¨qg©mª«g¬¤ ¢®j¯°l²±U³§ ´g¡gµx§ gug v¶lf¶hsik·ugv°_¸º¹ V ggWc§rgtjugv°l_fhi»k¼ugv»W l²½ ¾A©mª ¿c»pºÀ ´R¡RµdÁ c Âxªcgf+hji+kÄÃRsXYdjÅjÆRÇxw+y XRYxf+hji kÈuRv É ¯DlsÊAË \_^xw+y_ \_^xf+h i+k¼uv `RbRÌdc ÅjÆRÇxw+y_ `RbxV_Ógf+ÔWh Õgi+Ö k ugvUÃÍ ÎºÏxT c§´R ¡gµU ggRgugvAÐrgtUf+hºi¶kÑuRvqcgÒRcsRºdc Wc×WØj~gRÙ Í. 1. ,. 3. Nehari. (. 1).. ,. . Nehari. [8].. ,. ,. ,. ,. . PSfrag replacements θ. [5, 7].. S S1 Sn. Sn. θ S. PSfrag replacements. p1 θ. pn. θ. p2 S2 θ. PSfrag replacements θn S L1 r L 1 n Ln rn. r2 L2. θ2. D −→ PSfrag replacements L1 θ1. C1. C2. Cn. L2. edcrgtxf+hji+kÈuRvxw¶yÚRuRvWWcRgÛWj Vâ _Þgß ÜºÝ ö W÷cøé½gù8¾UyoË WªW gúgº xºcRàuá vdcº RØºRãjóDägá»²rtRåWuæ¶vdçéègcêUeRûqëlìx¢ íUdhmü îgjïRñuðxvxw°ògóxyÐrRcstô õu vWÞWcRUjTA rRtRuRvdw+yÚRuRvUWcTWRV Wjjý+©Úª þÀy§ÁRÿdcjàjá T ôRõWj ý l xp ql

(30) ºRuRvxw+yñrRtRuRvxWcRUjd cgº ©ÚÜºª Ý rgtRV uRvxw+yÚRugvWWcRWjdcgqloRRÛWjdcg§A¸ ¹»UjË y ÎÏUc§ g xcj u v»ºý©m ªÁ c eUc§r tUfh§ iÜjk¼ÝugvRc§ º â c§ôgõql Î àjá°l²½R¾D©mªW © ØºãgäR á°lmØxV TyÚT cj àá°lorgt f+hji+kÈuRvxw+V ymÎjÏdcgºRuRvWWcºRRÛWj» côRõ UØË U §»ú Ô g w ¿c ºcº Î

(31) rRtUf+hºi¶kÈugvxÃ µA f+hi+kÚl RüzR|gXR}gÂ "Î !WË U ÜjÍ Ý y_cºàjá V RWºdc °l_Î#dc$&%d'&(Uý Þ l) * ,+"-¤. /g¿cg 1: 3. .. ,. ,. ,. ,. [4].. ,. .. , Nehari [3]. ,. ,. .. ,. 1. 25. 3. ..

(32) ö á°©oªùDyñ â Øºãjägá°l76ØgË §Wp§Ájc T Í â Øãgägáq V ì R x c Uý Þ l²«g¬gv ¦»:9»ö g;8<=R §ãjäUÃ, =

(33) T 0"13 2845 ìc ö >g8? T ô õD©7/sã äxc@89q ¤Wl ©ÚªôRõ l &xË xp_àjá Í C"D&E Ãµqlo±3*gªRRö uRvdc¦RT¨ T ú V »¤s ÷ps, H8? 6&I: J&9¤§UT sãgädcRÞ Ád²ð»sògódc§ôRõ lK&L xj»cg& Mjª µdÃ *¶yÚw xý &N&O»,c PxWWH?sãä Ry QD©o, ª 9Dç§Ws: g×WË ©º Ã §ª rgtf¶hik¼ugvxc÷p « ¬gvg¦xcHSWÃ T&l7g U µW©_w¶ :H&V?Z V §dú T cjàjáqls¿cUúgú6WØË U VjT Îg#à VWYX : Vg3µ [Æqâ ©mª [_ÆUf¶hsiÞgßkml²n°p_ÎÏxc gu vUw¶y{ rgtxfhik·ugv c]gWjxcg§ : ©mª V ØãR ägáPgQg÷ SV gWdcºÜºàjÝ á°l 6 ØË *R V \8 Á cl T ØÍ Ë §àºUáqWl²÷»½ ¾ø ©m gªR RRuRvdcRRÛWjdccgàºáql â ØjãR÷ äRággWÛW÷ º öxáqcjà7l 6Wá»ØØUË UÃ T 0 132"4#5 " '8 A ( '&B '8( FG .. .. ,. 1. [6]. ,. ,. , .. ,. ,. ,. .. Joukowski. [1].. . . PSfrag replacements. ,. ,. ^`_badcfe. e N h gjilkbmondprqtsvuxw PRQgSxT ½R¾WË àjá TWV gRÛjUc"y#z{8|ÀºguRvWw+y{rRtguRvWc§Rjxcô õº xý Þ l) * ªM»çR}RÙUÃRÍ V T ÜjÝ c§½R¾UË sÞgß RÜºgÝ Ûjxc H&?dÁgÿ T Í ú~, *,ÞR/ß r tgugv»X YUf¶h§ik¼ug v+7l Ø°© H?°7l d y÷ gujvUw yÐrRtRuRvWc RWjdcºàjáql RË. 2. ,. .. , [4]. .. ,. ,. .. v. y. L1. L2. C1. C2. f (z). θ. θ. θ u. x Ln. Cn. RuRv XgYxf+hji+kÈuRv RuRvA²XRYxf+hºi+kÈuRv X}R~»c =[Æ l µDË ´¡µd RcRuv w+y / ÞXR}R~dc´R¡RV µd/ RRdcXRYdf+hji+ k¼uR¢v®g¯x8WM8cR Wjqlo Þ c§C « Þ8w+ y Ã W§cjUjpx8ý &8 cÁUÎ"»«¶ú CWc fhºi¶kÄc,@89AÐ£¤ V ºxý Þ §ÁR«¶ú T T Þ R XRYxl_f+Øxhjji+ª kÈuRvdcjÅjÆxf+hji+k¥Ã &AsgËql © l $&%d,'&(xý M»çRCWc & ºUý Þ V [Æ l²» XgYx f+hik Ë Uw+y >dc l ºË D. 2:. z. Jordan n. w = f (z). C1 , . . . , Cn. n. D. S. f (∞) = ∞, f 0 (∞) = 1, ,. D. S. Laurent. 0. (. , g(z), h(z). θ. , f (z). ,. f (z). C1 , · · · , Cn. 2).. D. f (z) = z + eiθ i(g(z) + ih(z)). ,. w. (1). L1 , · · · , L n. .. Im (e−iθ f ) = pm ,. z ∈ Cm , m = 1, · · · , n,. (2). 2. 26.

(34) &© V f+hji+k¥c,@&9qlxËRùxc« Þ T Í C»s¶l c V s C»¤xys ¢®R¯»&M& & && w+y/ , pm. ©mª "Å .ü. &. g(z). g(z) + Im(e−iθ z) = pm , z ∈ Cm , m = 1, · · · , n. ,. (3). g(z), h(z). ýql¤º »ü ¶y_#CjÛ/ xÃ&|dË xü Ë jxÀ dc ÃÞ 8 8& lºËx WÃjw C ö cºÎ"&OUw+y7/RWºdc V ºcW÷ps,$ %d, '&(xâ ý l)* +"-¤§ C > ØjãRäRáql_Øxºª c§ôRõqloYd»pC côgõ l g(∞) + ih(∞) = 0,. (4). g(∞) = h(∞) = 0. (4). (1). f (z). g(z), h(z). ,. g(z), h(z). f (z). F (z). F (z) = z + eiθ i(G(z) + iH(z)), Nl n X X. G(z) + iH(z) =. Qlj log (z − ζli ). MçRC» T & V cd# sãäR¯ T Í C"R« Þ M& gµ & l=1 j=1. , ζlj. Cl. Qlj. zmk. Nl n X X. Qlj log |zmk − ζli | − Pm = −Im(e−iθ zmk ),. ÚÎ"&Odc*Rc l=1 j=1. (5) (6). V R µ~» &¯. m = 1, 2, · · · , n, k = 1, 2, · · · , Nm. (7). & Nl X. Qlj = 0,. l = 1, 2, · · · , n. (8). ljËd÷ps«3* C rRtRugvdceRû» Íj".WÞª$&%d,'8(xý Þ c§«R¬ql¹»U T \_^x f¶hji+k·TuR v[ Z `Rbxf+hji kÈuRvjý¶©ÚªWÁRÿÎWc, Y8xw+y_ &|Wà&jql+"-xË UÃ j=1. ,. .. N gjilkbmoTWndV po`qtsvuxw e h ÜºÝ ÷ Þgß RgÛºdcºàºá zj|RXg}g~dc§r tguRvxc§¦»sã ä°l79 cH8?q ÿ:»V â Øjãgägá° lmØUjªggÛjxc§ôgT8õ°U lñYqp &§T Ã r tUf¶hºik¼ugvql²±U³ T Xg}g~ ãgä+l7J"9WË " * cugv ¦Uc V HSWP Ã ©m g ;x j\u ] v Ã g ÍUü rRtxf¶hji+kÈuRv»=T s nq psf+hi+k V c l_ØxR: ¡& ¢&Wjql W k_:£i+ k¥ÃgÔÎ"!xË V \&]ØË l_ØxUrRÃ txf+hji+ k¼uRvRÂ §Ô xl§Í rR\ t¤duRcv¦Rs¨RnAuRpsvWf+jhjRi+ W jxË ¥ W >dcx÷p§ Þ jxý l_Øx &© T &V f+hºi+k§¦g¯dcô õ Þ T V Ã c,& &V& :R×xË c ¢º®R¯ #c & V l jËW8Þ *Rcº« XYxf+hji+k uRvdc H&? T T s ôRõjWý l_ØUª ¨©Usð»sògó )* U_Ã. 3. ^`_badcfe. ,. , F (z). .. ,. .. 1. (. 3).. w. ,. .. ,. Joukowski. .. .. ψ(w) =. , W11 , W12. p p 1 (w + w − W11 w − W12 ) − ψ0 , 2 , ψ0. ψ(w). f (z). ,. . W11 , W12. w = F (z). ,. 3. 27. ,. , Joukowski. ,. ψ(w). D. (9). f −1. ψ0 = (−W11 − W12 )/4 [2]..

(35) rtxf°hji+V k]uRv c U÷ A lmuRv W© T&UÚl ±d³º zR|X}qV l RX }DsË V XR} ~ T c¦g¨sãRä°lJ&9xË cªÇd,H&SUÃ Ë «x ,c Rµ \ ¤ s c ¿cs¨»ÿsÊ\¤g~xÒ"¬xãRä ¯ l9YRü÷j +w y WcôRõRgWj l S. z. ∗. D. C1∗. D∗. ψ(w). , D∗. ∗. z∗ , D∗. .. ∗ ∗ ζ11 , . . . , ζ1N 1. ,. ∗. . D∗. .. D. ∗. z = F (z ). F ∗ (z ∗ ) = z ∗ +. ,Mç. . F ∗ (z ∗ ). . ψ(w). c,?&xý Þ l_Øxjª. N1 X j=1. , f −1 (w). loôRõxË . f −1 (w) ∼ Fe (w) = F ∗ (ψ(w)) = ψ(w) +. R« Þ &. V R µ. Q∗1 , . . . , Q∗N1. &. l j Ë d÷Wps«3* &© W÷ c. N1 X j=1. F (z1k ) PSfrag . replacements. p1. V. (11). j=1. ψ(w). lo ú T FGq©s ý Þ D. D. C1∗ ∗ ∗ ζ12 ζ11 ∗ ζ1N1. ψ(w). θ. ∗ Q∗j log(ψ(w) − ζ1j ).. (12) ψ ∗ (w). PSfrag replacements D∗. S. L1. ,. ∗ ∗ Q∗j log(z1k − ζ1j ) (k = 1, . . . , N1 ),. ∗ , z1k = ψ ∗ (F (z1k )). .. N1 X. (10). ,. ∗ ∗ z1k = Fe (z1k ) = z1k +. PSfrag replacements. ∗ Q∗j log(z ∗ − ζ1j ). z11. F ∗ (z ∗ ). C1. −→. −→. z12. z1N1. Wj l_Øx àºá rRtRuR\&vx] Ãd\_^xf+h i+k¼uv T Í H&? V c,&Al®¯Ë }RÙxÃÍ ¥ Ô W V eRûdc l²ugTvU´ c88O»sT×± °oª,²8ø³?#. c8°l ü exc§rRtUf¶hi¶k¼ugv V ÿÎcºà á pjWÃ TWV ÞV Þ ÎjjÏd c RcxR÷Rpsu,v H8?d rRtxf°cUhj÷Ri+pskÈ,uR¡8v¢8sWnAºps+f+7l hgWi+Øk¥Ë c àºá+uRlsv ¿RcWcúR R¶ú 6óØsË ×µ U°Úª zV · ©§ ¥ Ô x V ¶l xË Us · g \ ] l ¹2Rø ºdËdU T z Þ cRf Ã ¸ hºi¶kÚl»> ¼jªW»çWÁ T Ã f+hi¶k{c Þ Ã,½ R Úºxý Þ c&x c*R}gÙd c,§H& ?RV ¾Ç z jÃÞ cg´ f+¿hjsi+ z kmÀlo nq psRú uR vxf+w+h yÐrRi+txÁk ¦f+¯dhjci+ôkÈõduRcvWòWócRÁ Â=WjËdcgU»sÃ& Ã&MR<+jªW¤§Á $&ÄxË ¸ ·xT TUÍ V Ë T W¾X V Þ z Þ c [sÆxf+hjÞ i+kÚlonqpsRuRvql ´ p * /jugvRÂ Ô l »Ë T c °Rç /f°hji+k{c , × ©§§z cUgq l UØWË à áA lm½¾8©_ªW ¿W ºcjàjáqlñ RRRRuRvdcgRÛWjdcg6ØË Uglo »ú T ,c 3U Å ºuvql rRtdf°hºi+k uv §nA p P cgf°hji+kÚl f¶hºi¶Æk ¦R¯xc§ôgõql_¿j Çj ÚË Îcgf¶hi¶k¥c H&?Uc . 3:. ¡&¢&. ψ(w). , ψ(w). .. p p 1 (w ± w − W11 w − W12 ) − ψ0 2 . ψ(w) ,3. ψ± (w) =. ,2. (13). .. ,. .. D. ψ(w). .. ψ. . Joukowski. ,. ψ± .. .. .. [1].. ,. ,. .. ,. Ln ,. ,. ,. D, S W11 , W12 , . . . , Wn1 , Wn2. n. .. L1 , . . . , ψ. 4. 28.

(36) ÿ: P cgf¶hi¶k Þ l7»ËgjxcºjWý xË UÃ T W÷»ø cjj¶¤

(37) 8c [sÆD l XR}AË RRd c H&? ¡&¢&Wj l_Øx â ØjãRäRá Í V 1. . ψl. 1. z. l È ymw" [sÆ °© \8y²]x:& 8& l l:£»ik¥c,É8¦qÚË åWæjª, RP ¦xcgf+hji+k V c¦ ¨dc [sÆxf¶hji+kÐj¶¤§ c [sÆ»f+hji+kÚlon»»pouRvAl lo±x³z|RXR} lWØË àºá V XR}&MRºªzR|& Þ 01Y2&4"5 ì Cl∗. Ll (l ∈ {1, . . . , n}) ψl Wl1 , Wl2. (l). , Ll n−1. Cl. ψ(w). z∗. .. , Joukowski. . D. (l). D ,D. (l). XR}&MRjª jUý Þ l7WØq©sså&Êxý Þ _l Øx â ØjãRäRáA:ËgËxWÃ T RR dc,H&?Á +ÿ loå8Êxý Þ Ë ´ â ØjãgäRá° T lo xü z&Àd¡&¢8WºÞq0lÚ13Øx2& 4" 5U¶ç{z Þ f+hi+kÚl²nq p Rugvq7l Ì8ø pjU_Ã

(38) Uúø czg&| ìWc Í#?&.ÎÐÏºË ¨ Çql cå Êxý Þ cT Í"?&. T &xË x pjU T Í àºá » > , j / ¶ f h i k Á ¦ ¯ c 8 N x O c ° ² l C g g x c 8 H » ? : 8 ¡ ¢ º Ú l U Ø TWV f+hi+ké,c ¦ ¯ V W÷»ø\ ¤ ~d¤§ T ,c &N c , 8 j ¶ f º h i k Ñ ¦ ¯ V Ó °Ú Îgà gggxcjgugv»MgºªjØU å8ÊUý Þ V ¦g¯xc8N8O O»8Ò , cgT"

(39) jª V Ô Á Þ Õ&Ö ©_ ¿x ôRõWjxý &V NOql . *W cTWV Í?W" .ql_Ød C» c &N&Odc§ Í"?W" .d,c j W² Ã ×"ØÐjÃ /»R , ,. ∗ ∗ log(z ∗ − ζ11 ), . . . , log(z ∗ − ζ1N ) 1. ,w. (14). ψ. ∗ ∗ log(ψ(w) − ζ11 ), . . . , log(ψ(w) − ζ1N ) 1. .. (15). ,. (l). (l). log(ψl (w) − ζ1 ), . . . , log(ψl (w) − ζN ),. l = 1, . . . , n. (16). ,. .. , (10). (16). .. ψ(w). ,. C1∗. ψ(w) . ,. ∗. ∗. , F (z ). (16). .. ,. ψ1 (w), · · · , ψn (w). n. J(w) =. 1X ψl (w) n. (17). l_Øx WRË C"R~+÷»øÄôRõWjxý Þ lÙ&xË l=1. f −1 (w) ∼ Fe(w). = J(w) + eiθ i (G(z) + iH(z)). G(z) + iH(z) =. Nl n X X. Ã Ú y² öCÛ&g«& Þ V XRT}g~ d &g¯ |Wà"&dc ©Ú ª «3* U_Ã. /. (l). Qli log(ψ (l) (w) − ζli ). l=1 i=1. z (l). (l). (l). z11 , . . . , znNn (l = 1, . . . , n). .. (l). zmk = Fe(zmk ) = J(wmk ) +. Nl n X X. (l). (18). (l). Q∗lj log(zmk − ζj ),. (19). W÷»ø³>dc. m = 1, . . . , n, k = 1, . . . , Nm ,. (20). V l²8d ú T F Gq©ºý Þ g T ÷ c ú \_^/ `RbUf+ hºi¶kÈuRvWcsRj»"

(40) WºªWÁ ÜjÝ y_c§½R¾ ©Î «& W ÷dø_ÿ dcjàjáxÃ yo C (l). l=1 j=1. wmk = F (zmk ), zmk = ψl∗ (F (zmk )) F (zmk ) .. ψl (w). ,. 5. 29. D(l). ψl∗ (w).

(41) N\ ] §ådæ°çÐå Êdý Þ lØ» ª â Ø ãäáW÷ §Þß TÛWV »Þgcgß à áAl F=Gq© g Þ ß ggR ugvdc g ÛWºxcg»A©m ªWÁ:6Ø¶l¶×Ø»© à8á / &â 8ãR°l ©_gÃ+y{½R¾ ©sàjádc& ¡ ä8Oql å&æxË C ÜÞÝ ß. 4. Joukowski. ,. ,. .. çjèêéìë [1] Okano, D., Terazono, M., Ogata, H. and Amano, K.: Numerical conformal mapping from domains with multiple slits onto the canonical slit domains by the charge simulation method, Information, Vol.6, No.1, pp.107–118 (2003). [2] Okano, D., Ogata, H., and Amano, K.: Stagnation Point Analysis by the Numerical Conformal Mapping Using the Charge Simulation Method, Theoretical and Applied Mechanics – Proceedings of the 50th Japan National Congress on Theoretical and Applied Mechanics,. ÜºÝ Ù í Ý A îà#ï8ð X&ñ¾ò ógÉ88ô â Øãgägá»÷ ´g¡gµxs ggggu vxc Î ÞRß gWjdcjàjá õ8ö&¢& Ö"÷&øù"ú / / / Ü§Ý Ù â ØãºäºádmåWæ»ç³ûRà üg Þjß jRcà§á õ#ö"¢" Ö,÷#øùú ö ÓRÔÕRÖ / P &ø / ýþ ÿ zR| «gY â ØjãRäRáA_¿c×W&Ø /

(42) / / ýþ ÿ ÓRÔWÕRÖ ¡ / /. Vol.50, pp. 303–309 (2001). [3]. [4]. ,. ,. , ,. ,. :. Vol.42 No.3 pp.385–395 (2001).. :. ,. , Vol.31, No.5,. pp.623-632 (1990). [5] [6] [7]. :. 1989.. :. :. 1983.. (. ). 1973.. [8] Nehari, Z.: Conformal Mapping, McGraw-Hill, NewYork, 1952.. 6. 30.

(43) 領域分割による SOR 法の収束性に関する考察畑上到，坂谷金哉（熊本大学大学院自然科学研究科）. 1. 緒言近年の計算機の発達により，大規模数値シミュレーションが行われるようになってきている．特に，分散，並列型の計算機環境の登場は，さらに大規模な計算を可能にしていくものと思われる．しかしながら，それらの計算環境において，その性能を最大限いかし得るためには，プログラムを書き換える等，その環境に最適な形にかえていくことが必要となる．例えば，最も頻繁に使われる反復計算については，プログラムばかりか，数値スキームの構築自体に改良を加えなければならず，その困難さから，単調なスキームが選択されているのが現状である．このような中で，本研究では PVM(Parallel Virtual Machine) を用いた分散メモリ型の並列計算に着目し，この場合に問題となる数値スキームについて考察した．具体的には，２次元ポアソン方程式の境界値問題を有限差分法で離散化し，SOR(Successive Over Relaxation) 法を用いて数値的に近似解をもとめる問題について，いくつかのスキームを挙げて数値実験を行い，「分散メモリ型の並列計算の威力を最大限利用できるような数値スキームを構築する上で必要な情報はどのようなものであるか」について考察した．特に並列計算では避けることのできない計算領域の分割が，SOR 法の収束性に与える影響について，データ受け渡しの部分における加速係数に着目して調べた．さらに，２次元問題で得られた結果を３次元問題に適用し検証した．. 2. 取り上げた問題本研究では，以下のような２次元平面内の正方領域 Ω = {(x, y)|0 < x < 1, 0 < y < 1} におけるポアソン方程式のディレクレ型境界値問題を考える．   ∂2u ∂2u + 2 = 2x(x − 1) + 2y(y − 1) 2 ∂y  ∂x u(x, y) = 0. (x, y) ∈ Ω. (1). (x, y) ∈ ∂Ω.. この問題の厳密解は，. u(x, y) = x(x − 1) × y(y − 1). (2). であるが，正方領域 Ω を x 方向，y 方向それぞれに間隔 h の等間隔に分割し，(xi , yj ) = (ih, jh)，u(xi , yj ) =. ui,j として，式 (1) のポアソン方程式を以下のように離散化する． ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1 + = −fi,j . h2 h2. (3). ここで，−fi,j = 2xi (xi − 1) + 2yj (yj − 1) である．この離散化された式（連立一次方程式）を反復解法で解く場合，例えば SOR 法を適用すると， (k+1). ui,j. (k). = (1 − ω)uij +. ω 2 (k+1) (k+1) (k) (k) (h fij + ui−1,j + ui,j−1 + ui+1,j + ui,j+1 ) 4. (4). のようになる．ここで ω は加速係数，k はステップ数である．本研究では，この計算領域を図 1 のように分割した場合について，いくつかの反復解法について，その解法の性質や加速係数の依存性等について調べた．. 1. 31.