Book of Abstracts

EASIAM 2011
The 7th East Asia SIAM Conference&

RIMS Workshop on Methods in Industrial and Applied Mathematics
Waseda University, Kitakyushu Campus, June 27-29, 2011
http://www.oishi.info.waseda.ac.jp/~easiam2011/

The 7th East Asia SIAM Conference & RIMS Workshop on Methods in
Industrial and Applied Mathematics
EASIAM 2011

Kitakyushu Campus, Waseda University, Japan
June 27-29, 2011

Book of Abstracts

Organaizer

East Asia SIAM

Co-Sponsors
Research Institute for Science and Engineering, Waseda University
Research Institute for Mathematical Sciences (RIMS), Kyoto University
Numerical Analysis Group, Department of Computational Science andEngineering, Nagoya University
Kitakyushu City
West Japan Industry and Trade Convention Association

Scientific Committee
President
Hisashi Okamoto
Vice-President Ming-Chih Lai
Secretary
Shao Liang Zhang
Ex-President
Dongwoo Sheen(Kyoto University)(National Chiao Tung University)(Nagoya University)(Seoul National University)

Member(Universiti Brunei Darussalam)(Chinese Academy of Sciences)(Hong Kong Baptist University)(University of Hong Kong)(Bandung Institute of Technology)(Bandung Institute of Technology)(Ajou University)(Chonnam National University)(University of Macau)(University of Macau)(University of Malaya)(Universiti Sains Malaysia)(Ateneo de Manila University)(National University of Singapore)(National University of Singapore)(National Taiwan University)(Mahidol Univeristy)(Vietnam National University)(Hanoi University of Science)

Victor Didenko
Zhong-Ci Shi
Tao Tang
Wai-Ki Ching
L. H. Wiryanto
Hendra Gunawan
Hyung-Chun Lee
Byeong-Chun Shin
Xiao-Qing Jin
Haiwei Sun
Mohd Omar
Abdul Rahni
Rafael P. Saldana
Weizhu Bao
Zuowei Shen
I-Liang Chern
Yongwimon Lenbury
Pham Ky Anh
Nguyen Van Huu

1

GreetingThe 7th East Asia SIAM Conference & RIMS Workshop on Methods in
Industrial and Applied Mathematics will be held at Kitakyushu Campus,
Waseda University, Japan. The EASIAM conference gathers researchers inapplied mathematics with respect to science, engineering and technology.
The scope of the EASIAM conference covers all aspects of applied mathematics to science, industry, engineering and technology in the East and South
East Asia. Previous EASIAM meetings were held in Hong Kong (1st), Sapporo (2nd), Xiamen (3rd), Daejeon (4th), Brunei (5th) and Kuala Lumpur(6th). EASIAM is a section of Society for Industrial and Applied Mathematics (SIAM). It promotes basic research and education in mathematics, which
supports science and technology.
Topics of interest include applications of partial diﬀerential equations;applied and numerical linear algebra; ﬂuid mechanics; mathematical analysis of engineering problems; mathematical aspects of geophysical problems;
mathematical biology and/or medicine; mathematics of material sciences;
numerical analysis; operations research; optimization.
We are grateful to the distinguished SIAM keynote speakers & invited
speakers, who have accepted to give plenary talks:
SIAM Keynote Speakers• Mitchell Luskin (University of Minnesota, USA)
• Andrea Bertozzi (University of California, Los Angeles, USA)
Invited Speakers• Ji Hui (National University of Singapore, Shingapore)• Yasumasa Nishiura (Hokkaido University, Japan)• Tony Wen-Hann Sheu (National Taiwan University, Republic of China)• Roslinda Nazar (National University of Malaysia, Malaysia)• Eun-Jae Park (Yonsei University, Korea)
We wish you a productive, rich, instructive edition of the EASIAM conference and we thank you for your participation.

Hisashi Okamoto
President of the East Asia section of SIAM

Shin’ichi Oishi
Chair of local organizing committee
2

Acknowledfements
It is my great pleasure to thank all the participants of the 7th East
Asia SIAM conference, whose latest achievements are reported here. I am
particularly grateful to them since they do that when our country faces anunprecedented diﬃculty caused by the earthquake in March 11th, 2011, andthe nuclear plant disaster right after that.
On behalf of the committee of East Asia SIAM, I thank very much the
keynote lecturers and the invited speakers. Their determination to give a talkin Japan was vitally important in holding the conference in this hard time.
Last but not least, I would like to express my deep gratitude to organizers
of Waseda University, particularly Professor Shin ichi Oishi, without whose
eﬀorts this conference would not be held.

Hisashi Okamoto
President of the East Asia section of SIAM
Kyoto University

On behalf of the local organizing committee, we would like to express
out sincere thanks to Kitakyushu city for supporting the conference, and weare also grateful to Research Institute of Science and Engineering, Waseda
University for supporting the conference.

Shin’ichi Oishi
Chair of local organizing committee
Waseda University3Local Organizing Committee
Chair
MemberShin’ichi Oishi
Takeshi Ogita
Katsuhisa Ozaki
Xuefeng Liu
Naoya Yamanaka
Akitoshi Takayasu(Waseda University)(Tokyo Woman’s Christian University)(Shibaura Institute of Technology)(Waseda University)(Waseda University)(Waseda University)

4Free Shuttle Bus Information
Free shuttle buses are arranged for all conference participants.
Direction
Hotel (Comfort Hotel Kurosaki) ↔ Conference (Collaboration Center)
Date
27th, June
28th, June

29th, June

Departure TimeTime
Direction
7:45–8:00 am
Hotel→ Conference
19:45–20:00 pm Conference →Hotel
8:30–8:45 am
Hotel→ Conference
17:40–18:00 pm Conference → Banquet Hall
After Banquet Banquet Hall →Hotel
8:30–8:45 am
Hotel→ Conference
19:15–19:30 pm Conference →
Hotel

5

Public Transportation

No.
1
2
3
*4
5
6

Direction
Kitakyushu Airport ↔ Conference Place
Kitakyushu Airport ↔ Kokura Sta.
Kokura Sta.
↔ Kurosaki Sta.
Kurosaki Sta.
↔ Conference Place
Fukuoka Airport↔ Hakata Sta.
Hakata Sta.
↔ Orio Sta.

78

Orio Sta.
↔ Conference Place
Orio Sta.
↔ Kurosaki Sta.
* For more information, see next page.

Way
Kitakyushu Airport Access Bus
City Bus
Kagoshima Main Line
Public Bus
Fukuoka City Subway
Kagoshima Main Line
Limited Express Kirameki
Public Bus
Kagoshima Main Line

Taxi
Direction
Fukuoka Airport↔ Conference Place
Kitasyushu Airport ↔ Conference Place
Hakata Sta.
↔ Conference Place
Kokura Sta.
↔ Conference Place
Kurosaki Sta.
↔ Conference Place
Orio Sta.
↔ Conference Place

6

Fare (about)
13,000 JPY
14,000 JPY
19,000 JPY
8,000 JPY
4,000 JPY
2,000 JPYTime
80 min40 min20 min30 min5 min
47 min30 min20 min4 min

Fare
700 JPY
600 JPY
270 JPY
270 JPY
250 JPY
910 JPY
1910 JPY
200 JPY
200 JPY*4. Kurosaki Station ↔ Conference Place

Direction
Kurosaki Station → Gakkentoshi Hibikino (Conference Place)
Hour
8
9
10
11
12
13
14
15
16
17

Departure Time(Free Shuttle Bus Available)
47
12
52
22
52
07







Direction
Gakkentoshi Hibikino (Conference Place) → Kurosaki Station
Hour
9
10
11
12
13
14
15
16
17
18
19

Departure Time
48
12
47
17
02
18
18
(Free Shuttle Bus Available)





7

Banquet

Date
28th June, 19:00–21:30
Place
Hotel Harmonie Cinq ( TEL : 093-592-5401 )Address 12-3 Otemachi Kokurakita Kitakyushu Fukuoka, Japan, 803-0814.
• From Conference Place– Free Shuttle Bus (Direct only)
Departure Time : 18:00 on 28th June @ Collabolation center• From JR Kokura Station (2km away)
1. Bus No.27, No.110. Bus Stop: ”Soleil-hall-move mae”, 10min.
2. Taxi 10min, about 1,000 JPY.

8

10:40--

13:20--

15:20

15:20

15:40--

10:30

Lunch

1

2

3

General Session
29-3-D

General Session
29-2-D

Room D

Conference Center

General Session
29-3-C

General Session
29-2-C

Collaboration
Room C
Center

General Session
29-3-B

General Session
29-2-B

Keynote Lecture
29-1-A

13:20--

Room BMedia Center Room A

12:00--

29th, June
10:40--

Lunch

Conference Center

9:30--

15:20

4

General Session
29-4-D

General Session
29-4-C

General Session
29-4-B

15:40--

General Session
28-4-D

General Session
28-3-D

Room D

9:00--

General Session
28-4-C

General Session
28-3-C

General Session
28-4-B

Collaboration
Room C
CenterInvited Lecture28-2-A (1)(2)

15:40--

General Session
28-3-B

Invited Lecture
28-1-A

13:20--

Room BMedia Center Room A

12:00--

10:40--

28th, June

9:30--

Lunch

Conference Center

9:00--

General Session
27-4-D

Student Prize
27-3-A

Room D

12:00--

General Session
27-4-C

Invited Lecture27-2-A (1)(2)4Collaboration
Room C
Center

10:30

10:30

3

General Session
27-4-B

Opening

Media Center Room A

9:30--

Keynote Lecture27-1-A

2

Room B

9:00--

27th, June

1

Dinner

18:00--

BanquetTime

18:00--

Reception

18:00--

29th, June 9:00--10:40 @RoomA
10:40--18:00 @Conference Center

28th, June 9:00--12:00 @RoomA
12:00--18:00 @Conference Center

27th, June 8:30--15:40 @RoomA
15:40--18:00 @Conference Center

Registration Desk

Program

10

Monday, June 27, 2011
Opening Ceremony

09:00∼09:20

Room A

Chair: Shao-liang ZHANG◦ Opening talk by the president◦ Some remarks

27-1-A

SIAM Key Note Speaker

09:30∼10:30

Room A

Chair: Hisashi OKAMOTO◦ Mitchell Luskin : Atomistic-to-Continuum Coupling Methods . . . . . . . . . . . . . . . . . . . . . p.24

27-2-A(1)Invited Speaker

10:40∼11:20

Room A

Chair: I-Liang Chern◦ Wen-Hann Sheu : On a symplectic scheme that optimizes the dispersion-relation equationof the Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.26

27-2-A(2)Invited Speaker

11:20∼12:00

Room A

Chair: Xiao-qing Jin◦ HUI JI : Sparsity-based regularizations for image restoration . . . . . . . . . . . . . . . . . . . . . . p.2711

27-3-AStudent prize winner presentations

13:20∼15:20

Room A

Chair: Dongwoo Sheen◦ Yongyong Cai, Weizhu Bao and Hanquan Wang : Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates . . . . . . . . . . . . . p.28
◦ Lei Du, Tomohiro Sogabe and Shao-Liang Zhang : A variant of the IDR(s) method withquasi-minimal residual strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.30
◦ Chiun-Chang Lee : Limit Problems of Solutions for the Coupled Nonlinear Schrödinger
Equations and Steady-state Solutions of the Poisson-Nernst-Planck Systems . . . . . . . . . . p.32
◦ Akitoshi Takayasu and Shin’ichi Oishi : A method of computer assisted proof for nonlinear
two-point boundary value problems using higher order finite elements . . . . . . . . . . . . . . . . p.33

27-4-B

Numerical Linear Algebra

15:40∼17:40

Room B

Chair: T. Miyata and K. Ozaki◦ Takafumi Miyata, Lei Du, Tomohiro SOGABE, Yusaku YAMAMOTO and Shao-LiangZHANG : A Projection Approach to Complex Eigenvalues of a Specified Absolute Value
p.34◦ Takeshi Fukaya, Yusaku Yamamoto and Shao-Liang Zhang : Automatic Performance
Tuning for the Blocked Householder QR Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.36
◦ Xiao-Qing Jin, Wei-ping Shen and Chong Li : A Ulm Method for Inverse Eigenvalue
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.38
◦ Yung-Ta Li : Folding a linear system to a second-order form and its application to
structure preserving model order reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.39
◦ Roberd Saragih : Model reduction order of bilinear system using balanced singular perturbation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.41
◦ Katsuhisa Ozaki, Takeshi Ogita and Shin’ichi Oishi : Simplified semi-static floating-point
filter for 2D orientation Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.42

12

27-4-C

Mathematical Analysis

15:40∼17:40

Room C

Chair: C-A Wang and C-H. Cho◦ Ching-An Wang and Un-Un Kuo : On the existence of steady flows in a rectangularchannel with slip effect on one porous wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.44
◦ Yukihito Suzuki, Yukio Takemoto and Jiro Mizushima : Exchange of two modes of
instability leading to sustained oscillations in the flow past a cylinderical obstacle . . . . p.47
◦ Chien-Hong Cho and Marcus Wunsch : Global and singular solutions to the generalized
Proudman-Johnson equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.49
◦ Hendra Gunawan and Yoshihiro Sawano : Application of RKHS Theory to a Minimization
Problem with Prescribed Nodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.51
◦ Kwok Kin Wong, Robert Conte and Tuen-Wai Ng : Exact meromorphic solutions of thereal cubic Swift-Hohenberg equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.53

27-4-D

Mathematics for Flow

15:40∼17:40

Room D

Chair: L.H.Wiryanto and D. Victor◦ Leo Hari Wiryanto : Free-surface flow under a sluice gate of an inclined wall from deep
water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.55
◦ Syakila Binti Ahmad, Ioan Pop and Kuppalapalle Vajravelu : Free Convection Boundary
Layer Flow from a Vertical Cone in a Porous Medium Filled with a Nanofluid . . . . . . . . p.56
◦ Victor Didenko and Johan Helsing : Features of the Nystroem Method for the ShermanLauricella Equation on Piecewise Smooth Contours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.58
◦ Changtong Luo, Shao-Liang Zhang and Zonglin Jiang : A New Evolutionary Algorithm
for Ground to Flight Data Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.60
◦ Swaroop Nandan Bora and Santu Das : Ocean wave scattering within or over a porous
medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.62

13

Tuesday, June 28, 201128-1-A

Invited Speaker

09:30∼10:20

Room A

Chair: Hiroshi Kokubu◦ Yasumasa Nishiura, Takeshi Watanabe and Makoto Iima : Spatially localized travelingstructures and the asymptotic behavior in binary fluid convection . . . . . . . . . . . . . . . . . . . . p.64

28-2-A(1)Invited Speaker

10:40∼11:20

Room A

Chair: Hyung-Chun Lee◦ Eun-Jae Park : Adaptive mixed FEM and new Discontinuous Galerkin methods . . . p.66

28-2-A(2)Invited Speaker

11:20∼12:00

Room A

Chair: Mohd Omar
◦ Roslinda Mohd Nazar : Numerical investigation of mixed convection flow over a horizontal
circular cylinder and a solid sphere in nanofluids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.68

1428-3-B

Numerical Algorithm I

13:20∼15:20

Room B

Chair: H.Yasuda and C.Luo◦ Hidenori Yasuda : Lamellar pattern and two-phase shallow water equations . . . . . . . .p.70◦ Tanay Deshpande : Implementing the Finite Volume Method on a Rectangular Grid p.72
◦ Naoya Yamanaka, Masahide Kashiwagi and Shin’ichi Oishi : Accurate and Rigorous
Exponential Algorithm in Round to Nearest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.73

28-3-C

Phenomena, Model, Simulation I

13:20∼15:20

Room C

Chair: Tomoeda and C. Mendoza◦ Mohd Omar : An integrated manufacturing system for time-varying demand process p.75
◦ Takanori Katsura : Numerical analysis of a stationary transport equation by upwind
scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.76
◦ Akiyasu Tomoeda, Ryosuke Nishi, Kazumichi Ohtsuka and Katsuhiro Nishinari : Stable
Traffic Jam Occurs from Nonlinear Saturation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.78
◦ Carolina Mendoza and A.M. Mancho : Lagrangian Skeleton on the Gulf of Mexico Oil
Spill . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.80
◦ Keita Iida, Hiroyuki Kitahata, Satoshi Nakata, Masaharu Nagayama, : Mathematicalstudies on the self-motion of surfactant scrapings at the air-water interface . . . . . . . . . . . p.82
◦ Hisashi Okamoto : A curious bifurcation diagram arising in a scalor reaction-diffusion
equation with an integral constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .p.84

15

28-3-D

Applied Mathematics I

13:20∼15:20

Room D

Chair: K. Matsue and T. Mizuguchi◦ Hiroshi Kokubu : Application of topological computation method for global dynamics ofan associative memory model of the Hopfield type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.85
◦ Kaname Matsue : Rigorous verification of equilibria for evolutionary equations - Existence,
Uniqueness and Hyperbolicity - . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.86
◦ Tsuyoshi Mizuguchi, and Makoto Yomosa : Role of unstable symmetric solutions insymmetry restoring process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.88
◦ Edi Cahyono and Buyung Sarita : A Diffusion Equation Representing the Dynamics of
the Jakarta Composite Index (JSX) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.90
◦ Buyung Sarita, Femy Puspita Arisanty and Edi Cahyono : The Trend and Probability
Density Function of the Jakarta Composite Index (JSX) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.91
◦ Masataka Kuwamura : Dormancy of predators in prey-predator systems . . . . . . . . . . . p.92

28-4-B

Numerical Algorithm II

15:40∼17:40

Room B

Chair: M. Kimura and K. H. Kwon
◦ Masato Kimura and Kiyohiro Ishijima : Truncation error analysis for particle methods
p.94◦ Elliott Ginder, Seiro Omata and Karel Svadlenka : A numerical method for multiphasevolume-preserving mean curvature flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.96
◦ Min-Hung Chen and Rong-Jhao Wu : A High-Order Discontinuous Galerkin Method for
Elliptic Interface Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.98
◦ Teng-Yao Kuo, Hsin-Chu Chen and Tzyy-Leng Horng : A Coarse-Grain Parallel Scheme
for Solving Poisson Equation by Chebyshev Pseudospectral Method . . . . . . . . . . . . . . . . . . p.99
◦ Haiwei Sun and Qin Sheng : On the stability of an oscillation-free ADI method in laser
beam propagation computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.101

16

28-4-C

Phenomena, Model, Simulation II

15:40∼17:40

Room C

Chair: H. Suito and S.R. Pudjaprasetya◦ Hiroshi Suito, Takuya Ueda and Daniel Sze : Flows in thoracic aorta with torsion p.103◦ Alex Chang : A Pore-scale Network Flow Model for Two Phase Flow: The Genesis andMigration of Gas Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.105◦ Santiago Madruga and Santiago Madruga : Convective transport and stability in films ofbinary mixtures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.106◦ Sri Redjeki Pudjaprasetya, and Elis Khatizah : Longshore Wave Breaker with Reflected
Beach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.108◦ Ikha Magdalena and S.R. Pudjaprasetya : Wave Energy Dissipation in Porous Media
p.110◦ G Sharma and Madhu Jain : Multi-Compartment Modeling of Tumor Cells Interacting
With Dynamic Chemotherapeutic Drug . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.112

28-4-D

Applied Mathematics II & Networks

15:40∼17:40

Room D

Chair: K. Miura and W-K Ching◦ Matthew Min-Hsiung Lin : On LaSalle’s Invariance Principle and Its Application to
Synchronize Hyperchaotic Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.113◦ Keiji Miura : An unbiased estimator of noise correlations under signal drift . . . . . . . p.115◦ Wai-Ki CHING : Construction of Probabilistic Boolean Networks: A Maximum Entropy
Rate Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.117◦ Inseok Yang, Donggil Kim and Dongik Lee : Compensation Technique for Transmission
Delay and Packet Loss in Networked Control System based on Lagrange Interpolation p.119

17

Wednesday, June 29, 201129-1-A

SIAM Key Note Speaker

09:30∼10:30

Room A

Chair: Ming-Chih Lai◦ Andrea Bertozzi : Mathematics of Crime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.121

29-2-B

Numerical Algorithms III

10:40∼12:00

Room B

Chair: T. Ichinomiya and C-T. Wu◦ David Ni : Numerical Studies of Lorentz Transformation . . . . . . . . . . . . . . . . . . . . . . . . . p.122◦ Chin-Tien Wu, Ming-Chiang Jiang and Yu-Lin Tsai : Numerical Studies on Monge
Ampere Equation Arising from Free-Form Design of Geometric Optics . . . . . . . . . . . . . . p.124◦ Takashi Ichinomiya : Renormalization approach to solve Langevin’s equation . . . . . p.126◦ Kil H Kwon and Lee Jaekyu : Consistent approximation-sampling with multi pre andpost filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .p.128

29-2-C

Error Analysis

10:40∼12:00

Room C

Chair: K. Kobayashi and T. Kinoshita◦ Yintzer Shih and R. Bruce Kellogg : A Tailored Finite Point Method for Convection
Diffusion Reaction Problems with Variable Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .p.130◦ Kenta Kobayashi : On the interpolation constants over triangular elements . . . . . . . p.131◦ Takehiko Kinoshita, T. Kimura and M.T. Nakao : A posteriori estimates of inverse linear
ordinary differential operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.132◦ Xuefeng LIU , Akitoshi Takayasu and Shin’ichi Oishi : Numerical verification for solution
existence of elliptic PDE on arbitrary polygonal domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.133

18

29-2-D

Applied Mathematics III

10:40∼12:00

Room D

Chair: C. Y. Han and K.K. Viswanathan◦ Naoki Wada : Reconstruction of the density of tree-shaped networks from boundary
measurements of waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.135◦ Chang Yong Han and Song-Hwa Kwon : Cubic helical splines with Frenet-frame continuity
p.137◦ Viswanathan Kodakkal Kannan, Zainal Abdul Aziz and Saira Javed : Free vibrationof symmetric angle-ply laminated cylindrical shells of variable thickness including shear
deformation theory: spline method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .p.139

29-3-B

Bifurcation

13:20∼15:20

Room B

Chair: T. Sakajo and M. Shoji◦ Takashi Sakajo : Point vortex equilibria enhancing forces over two parallel plates . p.141◦ Tomoyuki Miyaji, Isamu Ohnishi and Yoshio Tsutsumi : Bifurcation analysis for the
Lugiato-Lefever equation in a square . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.143◦ Kuo-Chih Hung and Shin-Hwa Wang : A theorem on S-shaped bifurcation curve for apositone problem with convex-concave nonlinearity and its applications to the perturbed
Gelfand problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.145◦ Shin-Hwa Wang, Kuo-Chih Hung and Chien-Shang Yu : On the existence of a double
S-shaped bifurcation curve with six solutions for a combustion problem . . . . . . . . . . . . . p.146◦ Mayumi Shoji and Chika Shimizu : Appearance and disappearance of non-symmetricprogressive capillary-gravity waves of deep water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.148◦ Takeshi WATANABE, Makoto IIMA and Yasumasa NISHIURA : Bifurcation structureand spontaneous pattern formation in binary fluid convection . . . . . . . . . . . . . . . . . . . . . . . p.150

1929-3-C

Optimization I

13:20∼15:20

Room C

Chair: W.J.Leong and W-H Kuo◦ Wen-Hung Kuo : Single-machine group scheduling problems with time-dependent learning
effect and position-based setup time learning effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.152◦ Seyoung Park, Bong-Gyu Jang and Yuna Rhee : Optimal Retirement with Unemployment
Risks and Market Completion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.154◦ Wah June Leong, Mahboubeh Farid and Malik Abu Hassan : Rescaled Gradient-based
Methods with Fixed Stepsize for Large-scale Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . p.156

29-3-D

Numerical Analysis I

13:20∼15:20

Room D

Chair: T. Ishiwata and H. Notsu◦ Tetsuya Ishiwata : Behavior of polygonal curves by crystalline curvature flow with adriving force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.157◦ Guanyu Zhou and Norikazu Saito : Some remarks on the fictitious domain method withpenalty for elliptic problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.159◦ Norikazu Saito : Analysis of the finite volume approximation for a degenerate parabolic
equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.161◦ Takahito Kashiwabara : FEM analysis of the Stokes equations under boundary conditionsof friction type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.163◦ Hirofumi Notsu, Masahiro Yamaguchi and Daishin Ueyama : A mesh generator using aself-replicating system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.165

20

29-4-B

Accurate and High Performance Computing

15:40∼17:40

Room B

Chair: H. Fujiwara and D. Sheen◦ Hiroshi Fujiwara : A Remark on Numerical Instability of Complex Inverse Laplace Transforms using Multiple-Precision Arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.167◦ Dongwoo Sheen and Jiwoon Kim : Numerical Laplace inversion using mult-precision p.169◦ Chenhan Yu, Weichung Wang and Dan’l Pierce : CPU-GPU Hybrid Approaches inMultifrontal Methods for Large and Sparse Linear System . . . . . . . . . . . . . . . . . . . . . . . . . . p.170◦ Shugo Manabe : Design and implementation of a multiple-precision system on GPU p.172

29-4-C Optimization II

15:40∼17:40

Room C

Chair: G. M. Lee and M. Jain
◦ Gue Myung Lee : On Optimality Theorems for Robust Multiobjective Optimization
problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.174◦ Sangho Kum and Yongdo Lim : The Resolvent average on symmetric cones . . . . . . p.176◦ Dai-Ni Hsieh, Ray-Bing Chen, Ying Hung and Weichung Wang : Optimizing Latin Hypercube Designs by Particle Swarm with GPU Acceleration . . . . . . . . . . . . . . . . . . . . . . . . . p.178◦ Madhu Jain : Maximum Entropy Approach For Optimal Repairable Mx/G/1 Queue With
Bernoulli Feedback And Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.180◦ Afshin Ghanbarzadeh and Abbas : Application of the Bees Algorithm to Multi-Objective
Optimization Engineering Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .p.181

21

29-4-D

Numerical Analysis II

15:40∼17:40

Room D

Chair: C. Park and S-M. Chang◦ Chunjae Park and Dongwoo Sheen : A quadrilateral Morley element for biharmonicequations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.183◦ Shu-Ming Chang : Applying Snapback Repellers in Ecology . . . . . . . . . . . . . . . . . . . . . . p.185◦ Yusuke Morikura, Katsuhisa Ozaki and Shin’ichi Oishi : Verified solutions of linear systems
on GPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p.186

22

Abstracts

23

Atomistic-to-Continuum Coupling Methods
Mitchell Luskin
School of Mathematics, University of Minnesota
Minneapolis, MN 55455, USAe-mail: luskin@umn.edu http://www.umn.edu/~luskin/
Abstract
Materials often have highly singular strain fields at crack tips, dislocations, and grain boundaries that require the accuracy of atomistic modelingin small regions surrounding these defects. However, these localized defects
interact through long-ranged elastic fields with a much larger region where
the strain gradients are sufficiently small to allow accurate approximationby coarse-grained continuum finite element models. The size of the region
requiring atomistic modeling is typically several orders of magnitude smaller
than the size of the region that can be adequately approximated by coarsegrained continuum finite element models [6], which has motivated the development of numerical methods that couple atomistic regions with continuum
regions to compute up to length scales that are sufficiently large for accurate
and reliable scientific and engineering application.

Figure 1: The atomistic (red) region surrounding an edge dislocation is coupled to a continuum (grey) region. The O(1/r2 ) magnitude of the strain
gradient (where r is the distance to the dislocation) implies that a small
region (the core) around the dislocation requires atomistic modeling, buta coarse-grained continuum model can be accurately used outside the dislocation core. In the continuum region, the deformation of the atoms areapproximated by piecewise linear interpolation from the deformation of theblue atoms, thus reducing the computational degrees of freedom and allowing
significantly increased length scales to be computed.

24

During the past several years, we have developed an analysis of the stability, modeling error, and iterative solution of prototypical atomistic-tocontinuum coupling methods [1–3,5,6]. The formation and motion of defects
such as cracks, dislocations, or grain boundaries, occurs when the configuration loses stability, that is, when an eigenvalue of the Hessian of the energyfunctional becomes negative. It is thus essential for the evaluation of thepredictive capability of atomistic-to-continuum coupling methods that ourtheoretical analyses have investigated accuracy and stability for strains up to
the onset of instability of the atomistic energy [2, 6].
We have initiated a computational study of the accuracy of atomisticto-continuum coupling methods for several benchmark problems [4]. Our
theoretical analysis can successfully explain most of the computational results
for these multi-dimensional benchmark problems. However, we also observe
the need for additional theoretical analysis and benchmark problems to morethoroughly understand the accuracy of atomistic-to-continuum methods.
Joint work with M. Dobson, B. Van Koten, Xingjie Li, and C. Ortner.

References[1] M. Dobson, M. Luskin, and C. Ortner. Stability, instability, and errorof the force-based quasicontinuum approximation. Archive for Rational
Mechanics and Analysis, 197:179–202, 2010.
[2] M. Dobson, M. Luskin, and C. Ortner. Accuracy of quasicontinuum
approximations near instabilities. Journal of the Mechanics and Physics
of Solids, 58:1741–1757, 2010. arXiv:0905.2914v2.
[3] M. Dobson, M. Luskin, and C. Ortner. Iterative methods for the forcebased quasicontinuum approximation. Computer Methods in Applied
Mechanics and Engineering, to appear. arXiv:0910.2013v3.
[4] B. Van Koten, Xingjie Li, M. Luskin, and C. Ortner. A computational
and theoretical investigation of the accuracy of QC methods. In Num.
Anal. of Multi. Prob.. Springer, to appear. arXiv:1012.6031.
[5] M. Luskin and C. Ortner. Linear stationary iterative methods for theforce-based quasicontinuum approximation. In Numerical Analysis andMultiscale Computations, volume 82 of Lect. Notes Comput. Sci. Eng.
Springer Verlag, to appear. arXiv:1104.1774.
[6] B. Van Koten and M. Luskin. Development and Analysis of Blended
Quasicontinuum Approximations. arXiv:1008.2138.
25

On a symplectic scheme that optimizes the dispersionrelation equation of the Maxwell’s equations
Tony W. H. Sheu123*, L. Y. Liang1, C. M. Mei1, J. H. Li1
1. Department of Engineering Science and Ocean Engineering, National Taiwan
University, No. 1, Sec. 4, Roosevelt Road, Taipei, Taiwan, Republic of China
2. Taida Institute of Mathematical Science (TIMS), National Taiwan University
3. Center for Quantum Science and Engineering (CQSE), National Taiwan University

Abstract
In this talk I will present an explicit finite-difference time-domain scheme in Cartesian andcurvilinear coordinates for solving the Maxwell’s equations in non-staggered grids. The
proposed scheme for solving the Faraday’s and Ampere’s equations theoretically preserves thediscrete zero-divergence for the electrical and magnetic fields. The inherent local
conservation laws in Maxwell’s equations are also preserved discretely all the time using the
explicit second-order accurate symplectic partitioned Runge-Kutta temporal scheme. The
remaining spatial derivative terms in the semi-discretized Faraday’s and Ampere’s equations
are then discretized to provide an accurate functional relation that governs the numericalangular frequency and the wave numbers in two space dimensions. To achieve the goal ofgetting a better dispersive characteristics, a fourth-order accurate space centered scheme
which minimizes the difference between the exact and numerical dispersion relation equations
is proposed. Through the computational exercises, the proposed dual-preserving solver iscomputationally demonstrated to be efficient for use to predict the long-term accurate
Maxwell’s solutions.

26

Sparsity-based regularizations for image
restoration
Hui Ji
Department of Mathematics, National University of Singapore, 10 Lower
Kent Ridge Road, Singapore 119076
E-mail: matjh@nus.edu.sg
Abstract
Image restoration is about how to recover images distorted and degraded in
the formation or transmission stage. Most image restoration problems arechallenging ill-posed inverse problems that require certain prior informationof images to obtain satisfactory solutions. In recent years, the sparsity assumption of image/video data under suitable dictionaries, i.e. most of important information of images can be represented by only few coefficients in somesystem, has emerged as one powerful prior for most nature images. Combined with rapid progresses of `1 norm related minimization approaches forfinding the sparse approximation of a given signal in redundant systems (e.g.
split Bregman iteration), sparsity-based regularization has been one highly
efficient technique for solving many challenging image restoration problems.
In this talk, we will begin with a brief review on mathematical theory andnumerical algorithms of sparse approximation under redundant tight frame.
Then, several applications of sparsity-based regularization in image restoration will be presented, including: blind motion deblurring ([4]), blind image
inpainting ([1]), surface fitting in range imaging ([3]) and action analysis invideo ([2]).

References[1] B. Dong, H. Ji, J. Li and Z. Shen and Y.-H. Xu: Wavelet Frame Based Blind
Image Inpainting, Technical report, NUS, 2011.
[2] Y. Li, C. Fermuller, Y. Aloimonos and H. Ji: Learning shift-invariant sparse
representation of actions, IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), San Francisco, 2010.
[3] H. Ji, Z. Shen and Y.-H. Xu: Wavelet frame based scene reconstruction fromrange data, Journal of Computational Physics, 229 (6), March 2010, 2093–
2018.
[4] J. Cai, H. Ji, C. Liu and Z. Shen: Blind motion deblurring from a single image
using sparse approximation, IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), Miami, 2009.

27

Eﬃcient numerical methods for computing
ground states and dynamics of dipolar
Bose-Einstein condensates
Weizhu Bao1),2) , Yongyong Cai2),∗ , Hanquan Wang 4)1) Department of Mathematics, National University of Singapore, 117543,
Singapore
2) Center for Computational Science and Engineering, National Universityof Singapore, 117543, Singapore
3) School of Statistics and Mathematics, Yunnan University of Finance andEconomics, PR China∗ e-mail: caiyongyong@nus.edu.sg
Abstract
New eﬃcient and accurate numerical methods are proposed to compute
ground states and dynamics of dipolar Bose-Einstein condensates (BECs)
described by a three-dimensional (3D) Gross-Pitaevskii equation (GPE) witha dipolar interaction potential. Due to the high singularity in the dipolarinteraction potential, it brings signiﬁcant diﬃculties in mathematical analysis
and numerical simulations of dipolar BECs. In this paper, by decoupling
the two-body dipolar interaction potential into short-range (or local) andlong-range interactions (or repulsive and attractive interactions), the GPE
for dipolar BECs is reformulated as a Gross-Pitaevskii-Poisson type system.
Based on this new mathematical formulation, we prove rigorously existence
and uniqueness as well as nonexistence of the ground states, and discuss theexistence of global weak solution and ﬁnite time blowup of the dynamics in
diﬀerent parameter regimes of dipolar BECs. In addition, a backward Eulersine pseudospectral method is presented for computing the ground states
and a time-splitting sine pseudospectral method is proposed for computing
the dynamics of dipolar BECs. Due to the adaption of new mathematicalformulation, our new numerical methods avoid evaluating integrals with highsingularity and thus they are more eﬃcient and accurate than those numericalmethods currently used in the literatures for solving the problem. Extensive
numerical examples in 3D are reported to demonstrate the eﬃciency andaccuracy of our new numerical methods for computing the ground states anddynamics of dipolar BECs.

28

References[1] W. Bao, Y. Cai and H. Wang, Eﬃcient numerical methods for computing
ground states and dynamics of dipolar Bose-Einstein condensates, J.
Comput. Phys., 229 (2010), 7874-7892.

29

A variant of the IDR(s) method withquasi-minimal residual strategy
Lei Du1),∗ , Tomohiro Sogabe2) , Shao-Liang Zhang1)
1) Department of Computational Science and Engineering
Nagoya University
2) Graduate School of Information Science and Technology
Aichi Prefectural University∗ e-mail: lei-du@na.cse.nagoya-u.ac.jp
Abstract
We consider the solution of large, sparse and nonsymmetric linear systemsof the form
Ax = b(1)
where A is a nonsingular n × n real matrix, and b is a real vector of order n.
This problem arises in numerous applications such as discretizations of partial
differential equations, Newton-type methods for nonlinear linear systems andoptimization.
The IDR(s) method proposed by Sonneveld and van Gijzen [4] can solve
problem (1) effectively. Recently, it has received extensive attention and several variants have been developed. In this talk, we will describe a new variantof the IDR(s) method. The IDR(s) algorithm can converge faster than mostproduct-type methods, but also usually with irregular convergence behavior.
We reformulate the relations of residuals and their auxiliary vectors generated by the IDR(s) method in matrix form. Then, by this new formulationand motivated by methods in [1, 2, 3], we present the variant of IDR(s),named QMRIDR(s), to improve the demerit of its irregular convergence behavior. Both the fast and smooth convergence behaviors of the QMRIDR(s)method can be exhibited.
Two numerical results are reported to show the comparison between
IDR(s) and QMRIDR(s) in Figs. 1–2.

References[1] T.F. Chan, E. Gallopoulos, V. Simoncini, T. Szeto and C. Tong, A quasiminimal residual variant of the Bi-CGSTAB algorithm for nonsymmetricsystems, SIAM J. Sci. Comput., 15(1994), 338–347.

30

[2] R.W. Freund, QMR: a quasi-minimal residual method for non-Hermitian
linear systems, Numer. Math., 60(1991), 315–339.
[3] R.W. Freund, A transpose-free quasi-minimal residual algorithm for
non-Hermitian linear systems, SIAM J. Sci. Comput., 14(1993), 470–482.
[4] P. Sonneveld and M.B. van Gijzen, IDR(s): a family of simple and fastalgorithms for solving large nonsymmetric systems of linear equations,
SIAM J. Sci. Comput., 31(2008), 1035–1062.

Figure 1: PDE2961(n=2961, nnz=14585)

Figure 2: SHERMAN4(n=1104, nnz=3786)

31

LIMIT PROBLEMS OF SOLUTIONS FOR THE COUPLEDNONLINEAR SCHRÖDINGER EQUATIONS AND STEADY-STATE
SOLUTIONS OF THE POISSON-NERNST-PLANCK SYSTEMSCHIUN-CHANG LEE
DEPARTMENT OF MATHEMATICS, NATIONAL TAIWAN UNIVERSITY
Abstract. In this talk we present results from my Ph.D. thesis, which is joint works with myadvisor Professor Tai-Chia Lin and my Host Professor Chun Liu, when I visited the Institute ofMathematics and Its Applications (IMA) in Minnesota.
The Poisson-Boltzmann (PB) equation is conventionally used to model the equilibrium of bulk
ionic species in different media and solvents. We study a New Poisson-Boltzmann (PB n) equationin 1-dimensional interval (−1, 1):
2 φ00 =N2
X
bl βl e−bl φ
ak αk eak φ
−
R1
R1
a
φ
(y)
−bl φ (y) dykdy
l=1 −1 e
k=1 −1 e
N1Xwith a small dielectric parameter 2 and nonlocal nonlinearity which takes into consideration of
the preservation of the total amounts of each individual ion. This equation can be derived fromthe original Poisson-Nernst-Planck (PNP) system. Under Robin type boundary conditions withvarious coefficient scales, we demonstrate the asymptotic behaviors of one dimensional solutions of
PB n equations as the parameter  approaches zero. In particular, we show that in case of electroP N1
P N2
neutrality, i.e.,
k=1 ak αk =l=1 bl βl , we prove that φ ’s solutions of 1-D PB n equations may
tend to a nonzero constant c at every interior point as  goes to zero. The value c can be uniquelydetermined by ak , bl ’s valences of ions, αk , βl ’s total concentrations of ions and the limit of φ ’s at
the boundary x = ±1. Such a result can not be found in conventional 1-D Poisson-Boltzmann (PB)
P 1
P N2
equations. On the other hand, as N
k=1 ak αk 6=
l=1 bl βl (non-electroneutrality), solutions of 1-D
PB n equations may have blow-up behavior which also may not be obtained in 1-D PB equations.
We also study the compressible and incompressible limits of the solution of the non-linear
Schrödinger system. Recently, a rich variety of dynamical phenomena and a turbulent relaxation
have been observed in rotating Bose-Einstein condensates depicted by Gross-Pitaevskii equations
coupled with rotating fields and trap potentials. The dynamical phenomena range from shock-wave
formation to anisotropic sound propagation. The turbulent relaxation leads to the crystallizationof vortex lattices. To see the dynamical phenomena and the turbulent relaxation of two-component
rotating Bose-Einstein condensates, we study the incompressible and the compressible limits of twocomponent systems of Gross-Pitaevskii equations. We define ”H-function” a modulated energyfunctional which may control the propagation of densities and linear momentums under the effect ofrotating fields and trap potentials.

132

A method of computer assisted proof for
nonlinear two-point boundary value problems
using higher order finite elements
Akitoshi Takayasu1),∗ , Shin’ichi Oishi2,3)1) Graduate School of Fundamental Science and Engineering,
Waseda University2) Department of Applied Mathematics, Faculty of Science andEngineering, Waseda University
3) CREST, JST∗ e-mail: takitoshi@suou.waseda.jp
Abstract
Present authors have presented with Takayuki Kubo at University ofTsukuba a method of a computer assisted proof for the existence and uniqueness of solutions to two-point boundary value problems of nonlinear ordinary diﬀerential equations in the paper submitted for NOLTA, IEICE. Thismethod uses piecewise linear ﬁnite element base functions and sometimes
requires ﬁne mesh. To overcome this diﬃculty, in this paper, an improved
method is presented for the norm estimation of the residual to the operator
equation. In this reﬁned formulation, piecewise quadratic ﬁnite element basefunctions are used. A kind of the residual technique works sophisticatedly
well. It is stated that the estimation of the residual can be expected smaller
than that of the previous method. Finally, four examples are presented. Each
result demonstrates that a remarkable improvement is achieved in accuracyof the guaranteed error estimation.

References[1] S. Oishi: Numerical veriﬁcation of existence and inclusion of solutions
for nonlinear operator equations, Journal of Computational and Applied
Mathematics, 60 (1995), 171–185.
[2] M. Plum: Computer-assisted existence proofs for two-point boundary
value problems, Computing, 46 (1991), 19–34.
[3] N. Yamamoto and M.T. Nakao: Numerical veriﬁcations for solutions
to elliptic equations using residual equations with a higher order ﬁnite
element, Journal of Computational and Applied Mathematics, 60 (1995),
271–279.
33

A Projection Approach to Complex
Eigenvalues of a Specified Absolute Value
Takafumi MIYATA1),∗ , Lei DU1) , Tomohiro SOGABE2) ,Yusaku YAMAMOTO3) , Shao-Liang ZHANG1)
1) Graduate School of Engineering, Nagoya University
2) Graduate School of Information Science & Technology,
Aichi Prefectural University3) Graduate School of System Informatics, Kobe University∗ e-mail: miyata@na.cse.nagoya-u.ac.jp
Abstract
Given n × n matrices A and B, we consider computing a few eigenpairs(λ ∈ C, x ∈ C n ) satisfying
Ax = λB x( x ̸= 0 ).

This problem arises in many applications of scientiﬁc computing, such as
structural analysis and stability analysis of ﬂuid dynamics [1]. It is often
needed to ﬁnd a few eigenvalues near a speciﬁed point in the complex planeand for the problem several iterative methods are available [2, 3].
In recent years, so called the photonic crystals [4] have attracted much attention due to their potential applications in future integrated circuits basedon light, instead of electrons. In these practical problems, it is needed toﬁnd all the complex eigenvalues such that | λ | = 1, see Fig. 1. The numberof desired eigenvalues are small. However, they are usually scattered, notlocated near one point in the complex plane.
In this talk, we present an approach to complex eigenvalues of a speciﬁed
absolute value. This approach extends the projection method of Sakurai andSugiura [5] to a new method for the practical problems of photonic crystals.
In contrast to existing methods, our method can avoid computing almost all
the undesired eigenvalues and thus can be eﬃcient.

34

Im.

1

0

-1

-1

0

1

Re.

Figure 1: Distribution of eigenvalues in the complex plane. The desired
eigenvalues on the unit circle are shown by •, the other ones are shown by ◦.

References[1] Z. Bai, D. Day, J. W. Demmel, J. J. Dongarra, A test matrix collection
for non-Hermitian eigenvalue problems, 1996.
[2] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, H. van der Vorst, eds., Templates for the solution of Algebraic Eigenvalue Problems: A Practical
Guide, SIAM, Philadelphia, 2000.
[3] D. R. Fokkema, G. L. G. Sleijpen, H. A. van der Vorst, Jacobi-Davidson
style QR and QZ algorithms for the reduction of matrix pencils, SIAM
J. Sci. Comput., 20 (1998), 94–125.
[4] Y. Huang, Y. Y. Lu, S. Li, Analyzing photonic crystal waveguides by
Dirichlet-to-Neumann maps, J. Opt. Soc. Am. B, 24 (2007), 2860–2867.
[5] T. Sakurai and H. Sugiura, A projection method for generalized eigenvalue problems using numerical integration, J. Comput. Appl. Math.,
159 (2003), 119–128.

35

Automatic Performance Tuning
for the Blocked Householder QR Algorithm
Takeshi FUKAYA1),∗ , Yusaku YAMAMOTO2) , Shao-Liang ZHANG1)
1) Graduate School of Engineering, Nagoya University
2) Graduate School of System Informatics, Kobe University∗ e-mail: t-fukaya@na.cse.nagoya-u.ac.jp
Abstract
In developing software for scientific computations, it is important to optimize
the software in order to fully exploit the potential performance of hardware. Such
optimization has conventionally been done by humans based on their empiricalknowledge. However, as the computer architecture becomes complex, cost for optimizing becomes too high to be done only by humans. In addition, such optimization is very difficult for researchers in application fields who lack the knowledge
about hardware. These issues indicate the necessity of systematic and nonempirical frameworks for the performance optimization.
Under this situation, a framework named Software Automatic Tuning has recently been extensively researched [1]. In this framework, the target software isrequired to be tunable, that is, it has controllable variability or adaptability in it.
Additionally a tuning mechanism tunes the target software automatically so that itruns in good performance under given conditions. Considering these two points,
we have studied on the automatic tuning especially for the dense matrix computations such as LU or QR decomposition [2].
In this talk, we introduce our approach to automatic tuning for the blocked
Householder QR algorithm. The QR decomposition of A using the Householdertransformations is written in(I − tn y n y Tn ) · · · (I − t1 y 1 y T1 )A = R
where ti is a scalar, y i is a vector, and R is an upper triangular matrix. In theblocked algorithm, we first partition the target matrix into blocks and computepartial decomposition. Next, we aggregate the Householder transformations usedin the partial decomposition into a matrix form as(I − tp y p y Tp ) · · · (I − t1 y 1 y T1 ) = I − Y1 T1 Y1 T .
After that, we update the remaining blocks using matrix-matrix multiplication.
An overview of these computations is shown in Figure 1. Though we pay additional computation cost for aggregating the transformations, by using matrixmatrix multiplications, which have high data reusability, we can obtain advantage on recent architecture where floating point operation is much faster than datatransfer.
36






  

 







Figure 1: Overview of the blocked algorithm









Figure 2: Binary tree

The performance of the blocked algorithm depends on how to partition the target matrix into blocks. This is because total computation cost and the performanceof each matrix-matrix multiplication which affected by the matrix size are fixedby the partitioning way. Therefore, it is important to optimize how to partition thematrix. For this subject, optimization by man power has been done [3], however,there are few studies on systematic optimization.
Our goal is automatic tuning of the way of partitioning. In our approach,a binary tree shown in Figure 2 represents the way of partitioning. By using abinary tree as a parameter, various kinds of partitioning way can be controlled
systematically. Moreover taking advantage of the recursive structure in binary
trees, we construct the mechanism based on the dynamic programming to find the
optimal binary tree.
Through the performance evaluation on some computational environments, it
is shown that the performance obtained by our approach is as good as that obtainedby human tuning.

References[1] K. Naono, K. Teranishi, J. Cavazos, and R. Suda: Software Automatic Tuning: From Concepts to State-of-the-Art Results, Springer, 2010.
[2] T. Fukaya, Y. Yamamoto, and S. L. Zhang: A Dynamic Programming Approach to Optimizing the Blocking Strategy for the Householder QR Decomposition, Proceedings of IEEE Cluster 2008, 2008, 402–410.
[3] E. Elmroth, and F. Gustavson: Applying Recursion to Serial and Parallel
QR Factorization Leads to Better Performance, IBM Journal of Researchand Development, 44 (2000), 605–624.

37

A Ulm-like Method for Inverse Eigenvalue
Problems
Xiao-qing Jin1),∗ , Wei-ping Shen2) , Chong Li3) ,1) Department of Mathematics, University of Macau
2) Department of Mathematics, Zhejiang Normal University3) Department of Mathematics, Zhejiang University∗ e-mail: xqjin@umac.mo
Abstract
We propose a Ulm-like method for solving inverse eigenvalue problems,
which avoids solving approximate Jacobian equations comparing with otherknown methods. A convergence analysis of this method is provided andthe root quadratic convergence property is proved under the assumption of
the distinction of given eigenvalues. Numerical experiments are given incomparison with the inexact Newton-like method.
This paper is concerned with [1].

References[1] W. P. Shen, C. Li, X. Q. Jin: A Ulm-like method for inverse eigenvalue
problems, Applied Numerical Mathematics, 61(3) (2011), 356-367.

38

Folding a linear system to a second-orderform and its application to structure
preserving model order reduction
Yung-Ta Li
Department of Applied Mathematics, National Chiao Tung University,
Hsinchu, Taiwan 30010
e-mail: ytli@math.nctu.edu.tw
Abstract
A strategy that folds a linear system to a second-order representation ispresented in [1]. The strategy begins with the observation that the inverse
matrix of the transformation matrix permits a speciﬁc form. However, not
all matrices that promise the speciﬁc form are of full rank. Consequently one
might need to search for all possibilities. In other words, the constructionof the transformation matrix is combinatorially hard. We introduce a new
method to extract second-order system matrices by applying a similarity
transformation to a linear state space system. Instead of taking advantageof the speciﬁc form of the inverse of the transformation matrix, we obtainthe transformation matrix in a factored form. The new method solves thehard combinatorial problem for the construction of the transformation matrix
proposed in [1].
An application of the folding technique is for structure preserving reduction of large-scale second-order systems. The main idea is to reduce a linear
state space system and then fold it into a second-order form. A second-ordersystem with the output being a linear combination of displacement and velocity is considered in [1]. The output associated with displacement or velocity
only is considered in [2] or [3], respectively. Focusing on the folding process,we ﬁnd that the construction of the transformation matrices in [2] and [3]suﬀers the same challenge as in [1]. We will show that a combination usage ofthe new method and the procedure in [2] or [3] solves the hard combinatorial
problem.

References[1] D. G. Meyer and S. Srinivasan: Balancing and Model Reduction for
Second-Order Form Linear Systems, IEEE Trans. Aut. Control, 41(1996), 1632–1644.
39

[2] B. Salimbahrami and B. Lohmann: Order reduction of large scalesecond-order systems using Krylov subspace methods, Linear AlgebraAppl., 415 (2006), 385–405.
[3] B. Salimbahrami, B. Lohmann and A. Bunse-Gerstner: Passive reducedorder modelling of second-order systems, Mathematical and Computer
Modelling of Dynamical Systems, 14 (2008), 407–420

40Model Order Reduction for Bilinear System using
Balanced Singular Perturbation
Roberd Saragih
Industrial and Financial Mathematics Research Group, Faculty of Mathematics andNatural Sciences, Institut Teknologi Bandung, INDONESIA
e-mail: roberd@math.itb.ac.id
Abstract
Many problems in science and engineering have model in the bilinear systemsuch as electrical networks, mechanical links, nuclear fission, cardiovascularregulator, urban process, hormone regulation, predator-prey models, etc. In some
applications, it has the high–order degree of freedom (state). In such case, it isdifficult to implement the designed controller since it can cause the numericaldifficulties and high computational cost. Hence, it is desired the low-ordercontroller. Meanwhile, the low-order controller can be obtained from the low-ordermodel.
In this paper, we concerned with reduction order of the bilinear system modelby using the balanced singular perturbation approach. Some approximation
methods for the bilinear system, such as the balanced truncation and the projection
method, have been developed by authors [2,3]. The singular perturbation approachhas been applied to the linear system by Liu [1]. As given in [2], we first transformthe bilinear system into the balanced system. The concept of the singular
perturbation is based on the rejection of very fast modes from a model. Based on
this concept, we then divide the state of the balanced bilinear system into the fastmodes and the slow modes. The low-order model can be obtained by replacing thesmall parameter by zero to the fast mode. The advantage of the singular
perturbation is better than the balanced truncation in the low frequencies.

References
1]

Liu. Y. and Anderson, B.D.O. (1989), Singular perturbation approximation ofbalanced system, International Journal of Control, Vol.50, pp. 1379-1405.

2]

Liqian Zhang and James Lam, (2002), On H2 model reduction of bilinear
systems, Automatica, 38, pp. 205-216.

[3] Zhaojun Bai and Daniel Skoogh (2006), A projection method for model
reduction of bilinear dynamical systems, Linear Algebra and Its
Applications, 415, pp. 406-425

41

Simplified semi-static floating-point filter for
2D orientation problem
Katsuhisa Ozaki 1),4)∗ , Takeshi Ogita 2).4) , Shin’ichi Oishi 3),4) ,1) Department of Mathematical Sciences, Shibaura Institute of Technology2) Department of Mathematical Sciences,
Tokyo Woman’s Christian University
3) Faculty of Science and Engineering, Waseda University
4) Japan Science and Technology Agency, CREST∗ e-mail: ozaki@sic.shibaura-it.ac.jp
Abstract
This talk is concerned with robustness problems in computational geometry, especially 2D orientation problem is focused on. The 2D orientation
problem is one of basic problems in the ﬁeld of computational geometry. The
following is brief explanation of this problem: An oriented line and a point
are given on the two-dimensional surface. The aim is to clarify that the pointis left or right, or on the line. It can be boiled down to a sign of the 3-by-3
matrix determinant. If the sign of the determinant is positive / negative, thepoint lies on the left / right to the line. If the determinant is zero, the pointis on the line. Based on this problem, segment intersection and a point-in
polygon problem can be solved. Although the computation is small-scale, the2D orientation problem plays important roles in computational geometry.
Floating-point arithmetic is fast performed by modern architectures. However, the number of signiﬁcand bits of a ﬂoating-point number is ﬁnite so thatit may cause rounding errors in each operation. In the worst case, unreliable
results may be obtained due to accumulation of rounding errors. As for the
2D orientation problem, the computed determinant may be opposite in sign
to the exact determinant. The accuracy is especially worried when the pointis very close to the line.
Classical algorithms in computational geometry correctly work if the allcomputations are exact. Namely, it is assumed that all expressions are evaluated by rational arithmetic (symbolic computations). Once an inexact result
is obtained, then an algorithm may output a meaningless result. For examples, there are several algorithms outputting a convex hull which is theminimum convex enclosing all points. However, meaningless results are sometimes output, for example, a resultant convex does not contain all points or
the result is not convex. Such problems are called robustness problems andmany examples can be shown in [1].
42

If the accuracy of the computed result is a concern, using multi-precision
arithmetic or rational arithmetic is an option. However, these approaches
take much computing time. In addition, many problems are correctly solved
by ﬂoating-point arithmetic. Therefore, straightforward use of robust computations is not eﬀective.
There are so-called ﬂoating-point ﬁlters. They quickly guarantee that thecomputed sign of the determinant is correct if the problem is well-conditioned.
It is a kind of self-validating methods. If the ﬁlter cannot guarantee thecorrectness of the result, then robust computations should be applied. As
for earlier researches for the ﬂoating-point ﬁlters, for example, semi-static
ﬁlters [2, 3] and dynamic ﬁlters [4, 5] have been developed.
We suggest a simpliﬁed semi-static ﬂoating-point ﬁlter for the 2D orientation problem. When the ﬁlter handles ﬂoating-point exceptions like overﬂow
and underﬂow, then the checks of the order of the magnitude for input coordinates seem to be necessary. Several branches are required so that it makes
the performance slow down. Our ﬁlter handles any ﬂoating-point exceptions
with only a branch so that it is possible to work fast. Finally, numericalexamples are shown to illustrate the eﬃciency of the proposed ﬁlter.

References[1] L. Kettner, K. Mehlhorn, S. Pion, S. Schirra, C. Yap: Classroom Examples of Robustness Problems in Geometric Computations, Computational Geometry, 40, 61–78, 2008.
[2] J. R. Shewchuk: Adaptive Precision Floating-point Arithmetic and Fast
Robust Geometric Predicates, Discrete & Computational Geometry, 18,
305–363, 1997.
[3] G. Melquiond, S. Pion: Formally certiﬁed ﬂoating-point ﬁlters for homogenous geometric predicates, Theoretical Informatics and Applications, Special issue on Real Numbers, 41:57–69, 2007.
[4] V. Y. Pan, Y. Yu: Certiﬁed Computation of the Sign of a Matrix Determinant, Proc. 10th Annual ACM-SIAM Symposium on Discrete Algorithms, ACM Press, New York, and SIAM Publications, Philadelphia,715-724, 1999.
[5] H. Brönnimann, C. Burnikel, S. Pion: Interval Arithmetic Yields Efﬁcient Dynamic Filters for Computational Geometry, Discrete Applied
Mathematics, 109:25–47, 2001.

43

On the existence of steady flows in arectangular channel with slip effect on one
porous wall
Un-Un Kuo and Ching-An Wang∗
Department of Mathematics, National Chung Cheng University, TAIWAN∗ e-mail: cawang@math.ccu.edu.tw
Abstract
We study a boundary layer problem arising from the study of steadylaminar flows in channels with one porous wall with the slip effect. By
applying a similarity transformation [1], the governing equation of the NavierStokes system for steady, incompressible laminar flows in this channel could
be reduced to
f 000 (η) + R((f 0 (η))2 − f (η)f 00 (η)) = K.
(1)
Here K is an integration constant and R corresponds to the cross flow
Reynolds number based on wall velocity (filtration Reynolds number), while
positive (negative) R represents the suction (injection) through the walls.
The function f is related to the stream function, η is the normalized transverse coordinate, namely, η = 1 at the the porous wall, while the position at
the insulated wall is denoted by η = −1.
By imposing the slip effect on the porous wall, the steady flows could
be studied from the boundary value problem (BVP) of (1), subjects to theboundary conditions:f (−1) = f 00 (−1) = 0,

f (1) = 1,

f 0 (1) + ϕf 00 (1) = 0,(2)where positive ϕ is the slip coefficient.
Following the frame work in [9], it is our purpose to explore the solutionsof (BVP). It is clear that the given problem possesses the unique solution
when R = 0. For nonzero R, as in [2, 4], suppose we set f (η) = bG(ξ)/2R,
ξ = b(η + 1)/2 for some positive b which is to be determined. Then, G(ξ)
satisfies the following associated boundary value problem(BVP1 ):G000 (ξ) + (G0 (ξ))2 − G(ξ)G00 (ξ) = 16RK/b4 ≡ β,(3)G(0) = G0 (0) = 0,

(4)

G(b) = 2R/b,

G0 (b) + ϕbG00 (b)/2 = 0.

44

(5)

By assigning values to
G00 (0) = α,G000 (0) = β,(6)we let G(ξ; α, β) be the solution of the initial value problem (IVP) of (??),(??) and (??). Suppose, for a prescribed positive ϕ, G0 (b; α, β)+ϕbG00 (b; α, β)/2 =0 holds at ξ = b∗ , then (BVP) will possess a solution with R = b∗ G(b∗ ; α, β)/2
and K = (b∗ )4 β/16R.
From the homogeneity of G(ξ; α, β),R(α, β) = R(α/λ3 , β/λ4 ),(7)K(α, β) = K(α/λ3 , β/λ4 ),(8)for all λ > 0, we are able to verify that given any position ϕ > 0, (BVP) canonly possess either nonnegative and monotone, or non-monotone solutions.
In particular, there exists a positive number R(ϕ) such that (BVP) possesses
at least one nonnegative and monotone solutions for R < R(ϕ), and at least
non-monotone solution for R ≥ R(ϕ). Moreover, the given problem possesses
no solution for R ≤ 0, K ≥ 0.

References[1] A.S. Berman, Laminar flow in channels with porous walls, J. Appl.
Phy., 1953, 24, 1232-1235.
[2] F. M. Skalak and C. Y. Wang, On the nonunique solutions of laminarflow through a porous tube or channel SIAM J. Appl. Math., 1978, 34,
535-544.
[3] S. M. Cox, Analysis of steady flow in channel with porous walls, or withaccelerating walls, SIAM J. Appl. Math., 1991, 51, 429-438.
[4] T.-W. Hwang and C.-A. Wang, On multiple solutions for Berman’s problem, Proc. Roy. Soc. Edinburgh, 1992, 121A: 219-230.
[5] S. Chellaml, M. Wiesnerl and C. Dawson2,Slip at a unifromly porous
boundary: effect on fluid flow and mass transfer, Journal of Engineering
Mathematics, 1992, 26, 481-492.
[6] C.-A. Wang and Y.-Y. Chen, On existence of similarity solutions for
Laminar flow in a channel with one porous wall, Bulletin Inst. Math.
Academia Sinica, 1992, 20, 299-320.

45

[7] C.-A. Wang and T.-C. Wu, Similarity solution of steady flows in a channel with accelerating walls, Computers Math. Applic., 1995, 30, 1-16.
[8] S. Chellam and M. Liu, Effect of slip on existence, uniqueness,and behavior of similarity solutions for steady incompressible laminar flow in porous tubes and channels, Phys. Fluid, 18(8), 2006,
DOI:10.1063/1.2236302.
[9] U.-U. Kuo and C.-A. Wang, Multiple solutions of steady flows in a rectangular channel with slip effect on two equally porous walls, to appear
in Taiwanese J. Math.

46

Exchange of two modes of instability leading
to sustained oscillations in the ﬂow past acylinderical obstacle
Yukihito Suzuki1),∗ , Yukio Takemoto2) , Jiro Mizushima2)
1) School of Fundamental Science and Engineering, Waseda University
2) Faculty of Science and Engineering, Doshisha University∗ e-mail: y_suzuki@aoni.waseda.jp
Abstract
Flows past cylindrical obstacles often exhibit staggered arrangement oflocalized vortices, called Bénard-Kármán’s vortex street. The most typicaland simplest example is the ﬂow around a circular cylinder placed in a uniform ﬂow. The vortex street was conjectured to arise due to instability of theotherwise steady ﬂow and the classical stability theory, based on a parallel
ﬂow approximation, was applied to the ﬂow[1], which predicted a so smallvalue such as Rec = 3.2 for the critical Reynolds number of instability. Herethe Reynolds number is deﬁned with the cylinder diameter and the uniform
velocity. On the other hand, experiments and numerical simulations showthat the ﬂow becomes oscillatory for Re & 40 − 50. In fact, the global stability analysis of the ﬂow without resorting the parallel ﬂow approximation
revealed that the critical Reynolds number Reg is 46.184[2]. The diﬀerence
between Rec = 3.2 and Reg = 46.184 is beyond allowance of error arisingfrom the parallel ﬂow approximation in the classical stability theory.
In order to resolve the disagreement between the classical theory and theglobal stability analysis, a packet of disturbance is considered instead of a
monochromatic wave[3, 4, 5]. For disturbances having a packet shape, twodistinct instability modes are identiﬁed, one of which is called convective
instability, and the other absolute instability. The ﬂow is said convectively
unstable when a packet of growing disturbances has a positive group velocity,and the absolute instability is deﬁned if there exists a packet of growing disturbance with zero group velocity. Many numerical works have been reported
on the convective and absolute instabilities of the cylinder wake based on theparallel ﬂow approximation and it was found that the value Rec = 3.2 fromthe classical stability theory indicates the critical condition for the convective
instability, while the criticality for the absolute instability is determined asRea ∼ 25[4]. It is a common understanding at present that the global instability occurs after the region of the absolute instability has a greater extent
in the streamwise direction than 3.5 times the cylinder diameter[4] so that
there exists a ‘pocket’ between Rea ∼ 25 and Reg = 46.184[5].
47We aim to resolve the gap between the absolute and global instabilities
and clarify the mechanism of sustained oscillations in the ﬂow, in which impulsive force is added at a certain point in the ﬂow and subsequent spatiotemporal development of induced disturbance is observed. In order to clarify theunderlying physics of the transition from globally stable to unstable ﬂow, we
introduce the notions of passive and active modes of instability, extending
the convective and absolute instabilities to non-parallel ﬂows, respectively.
Our numerical analysis brought us a conclusion that the oscillation is sustained by superiority of the growth of disturbance due to instability over thedecrease due to advection of the packet. According to our deﬁnitions of thepassive and active instabilities, the ﬂow becomes unstable to active mode of
disturbance simultaneously in the whole ﬂow ﬁeld at the critical Reynolds
number Reg ∼ 47 for the global instability, whereas some extent of the regionis passively unstable above a critical Reynolds number Rep . When the ﬂow
is passively unstable, the ﬂow can be oscillatory if some external disturbanceis continuously added anywhere in the ﬂow ﬁeld. Otherwise, any localized
disturbance is swept out to far downstream at Re < 47. At supercritical
Reynolds numbers, any packet of disturbance added in any location will expand over the entire region behind the cylinder, proceeding upstream to theregion just behind the cylinder as well as being advected far downstream.
Once the packet of disturbance arrives at the vicinity of the cylinder, it activates oscillation leading to the global instability and becomes an oscillation
source.
We will propose an one-dimensional model equation, like a GinzburgLandau equation, which illustrates the mechanism of sustained oscillationand the transition from the passively unstable to actively (globally) unstable
ﬂow. The equation and its solutions will be presented at the conference.

References[1] S. Taneda, J. Phys. Soc. Japan., 18, (1963), 288.
[2] C. P. Jackson, J. Fluid Mech., 182, (1987), 23.
[3] H. Oertel, Jr, Annu. Rev. Fluid Mech., 22, (1990), 539.
[4] P. Huerre and P. A. Monkewitz, Annu. Rev. Fluid Mech., 22, (1990),
473.
[5] J. M. Chomaz, Annu. Rev. Fluid Mech., 37, (2005), 357.
[6] Y. Takemoto and J. Mizushima, Phys. Rev. E, 82, (2010), 056316.
48

Global and singular solutions to thegeneralized Proudman-Johnson equation
Chien-Hong Cho1),∗ , Marcus Wunsch2)
1) Department of Mathematics, National Chung Cheng University
2) Research Institute for Mathematical Sciences, Kyoto University∗ e-mail: chcho@math.ccu.edu.tw
Abstract
We consider the Generalized Proudman-Johnson equation

ftxx + f fxxx = afx fxx , x ∈ R, t > 0,
f (0, x) = f0 (x)(1)

on the whole real line as well as the periodic Generalized Proudman-Johnson
equation
ftxx + f fxxx = afx fxx , x ∈ S = R/Z, t > 0.
(2)
f (0, x) = f0 (x), x ∈ S
The generalized Proudman-Johnson equation (1), (2) can be motivatedin several ways. The first motivation comes from the study of the twodimensional incompressible Euler equations. For parameters a ∈ R ∪ {∞},
Eq. (1) was first derived by Okamoto and Zhu [7]. Remarkably, the parameter a interconnects several well-studied equations within the framework of
the generalized Proudman-Johnson equation (1), (2). For a = 1, it reduces
to the original Proudman-Johnson equation [8], which is obtained via theseparation of spatial variables ψ(t, x, y) = yf (t, x) for the stream function ψof the two-dimensional incompressible Euler equations, where x is containedin a finite interval, and y ∈ R. Moreover, the generalized Proudman-Johnson
equation encompasses the well-known Burgers equation for a = −3, which isa successful model in gas dynamics; if a = −2, the generalized ProudmanJohnson equation becomes the Hunter-Saxton equation [5] arising in the
study of nematic liquid crystals.
The existence of global and singular solutions were widely investigatedrecently, see for instance [1], [2], [6], [9]. In this talk, we construct a weaksolution for (1) for −2 ≤ a < −1 by the modified characteristic method.
Then by using an abstract lemma in [3], we explore the blow-up solutions for(2) for a < −1 and their exact blow-up rate.
49

References[1] A. Bressan, A. Constantin: Global solutions of the HunterVSaxton equation, SIAM J. Math. Anal., 37 (2005) 996-1026.
[2] X.-F. Chen and H. Okamoto: Global existence of solutions to the generalized Proudman-Johnson equation, Proc. Japan Acad. Ser. A, 78 (2002)
136-139.
[3] A. Constantin and J. Escher: Wave breaking for nonlinear nonlocal
shallow water equations, Acta Math., 181 (1998) 229-243.
[4] C.-H. Cho and M. Wunsch: Global and singular solutions to the generalized Proudman-Johnson equations, J. Differential Equations, 249 (2010)392-413.
[5] J.K. Hunter, R. Saxton: Dynamics of director fields, SIAM J. Appl.
Math., 51 (1991) 1498-1521.
[6] H. Okamoto: Well-posedness of the generalized ProudmanVJohnson
equation without viscosity, J. Math. Fluid Mech., 11 (2009) 46-59.
[7] H. Okamoto, J. Zhu: Some similarity solutions of the NavierVStokes
equations and related topics, in: Proceedings of 1999 International Conference on Nonlinear Analysis (Taipei), Taiwanese J. Math., 4 (2000)
65-103.
[8] I. Proudman, K. Johnson: Boundary-layer growth near a rear stagnation
point, J. Fluid Mech., 12 (1962) 161-168.
[9] R. Saxton, F. Tiglay: Global existence of some infinite energy solutions
for a perfect incompressible fluid, SIAM J. Math. Anal., 4 (2008) 14991515.

50

Application of RKHS Theory to a
Minimization Problem with Prescribed Nodes
Hendra GUNAWAN1),∗ , Yoshihiro SAWANO2)
1) Department of Mathematics, Institut Teknologi Bandung, Indonesia2) Department of Mathematics, Kyoto University, Japan∗ e-mail: hgunawan@math.itb.ac.id
Abstract
Let d2 < α < ∞. We define a Hilbert space Hα to be the set of functions
f on [0, 1]d of the formf (x1 , . . . , xn ) :=X

am1 ···md sin(m1 πx1 ) · · · sin(md πxd )m1 ,...,md ∈N

for whichkf kHα :=X
π 2α(m2 + · · · + m2d )α |am1 ···md |2 < ∞.
2d m1 ,...,md ∈N 1

The above norm is induced from the inner producthf, giHα =X
π 2α(m2 + · · · + m2d )α am1 ···md bm1 ···md ,d2 m1 ,...,md ∈N 1

where am1 ···md and bm1 ···md are the coefficients of f and g, respectively. We
are then interested in studying the following problem:
Minimize kf kHα
subject to the prescribed nodes:
f (pk ) = ck ,k = 1, . . . , N,

where pk := (pk1 , . . . , pkd ) ∈ (0, 1)d and ck ∈ R are given. [Note here thatthe points pk ’s are ‘inside’ the unit cube [0, 1]d .]The 1- and 2-dimensional case have been studied by Gunawan et al. [1, 2].
In this talk, we shall present the use the theory of reproducing kernel Hilbert
spaces to study the problem (in a more general setting). One of our results
is the following:
Proposition. The solution to the minimization problem
Minimize kf kHα51

subject tof (p1 , . . . , pd ) = 1,
is given by F (x1 , . . . , xd ) :=AX

sin(m1 πp1 ) · · · sin(md πpd )
sin(m1 πx1 ) · · · sin(md πxd ).
(m21 + · · · + m2d )αm1 ,...,md ∈N

where A−1 :=P
m1 ,...,md ∈N(sin(m1 πp1 )··· sin(md πpd ))2.
(m21 +···+m2d )αA more general result will be presented. In addition, we shall also consider
the Hölder continuity of the solution, by using properties of Besov spaces.

References[1] H. Gunawan, F. Pranolo and E. Rusyaman, “An interpolation methodthat minimizes an energy integral of fractional order”, in D. Kapur(ed.), ASCM 2007, LNAI 5081 (2008), 151–162.
[2] H. Gunawan, E. Rusyaman and L. Ambarwati, “Surfaces with prescribed nodes and minimum energy integral of fractional order”, submitted.
[3] G.G. Lorentz, Approximation of Functions, AMS Chelsea Publishing,
Providence, 1966.
[4] M. Sugimoto and N. Tomita, “The dilation property of modulation
spaces and their inclusion relation with Besov spaces”, J. Funct. Anal.
248 (2007), 79–106.
[5] M. H. Taibleson, “On the theory of Lipschitz spaces of distributions on
Euclidean n-space. I. Principal properties”, J. Math. Mech. 13 (1964),
407–479.

52

Exact meromorphic solutions of thereal cubic Swift-Hohenberg equation
Robert Conte1,2∗ , Tuen-Wai Ng1† and Kwok-Kin Wong1‡

1. Department of Mathematics, The University of Hong Kong.
2. LRC MESO, Centre de mathématiques et de leurs applications (UMR8536) et CEA-DAM, École normale supérieure de Cachan.
E-mail: Robert.Conte@cea.fr, ntw@maths.hku.hk, wkkm@graduate.hku.hk

Abstract
We considered the real cubic Swift-Hohenberg (RCSH) equation

2
∂u
∂2
∂2
= εu − 1 + 2 + 2 u − u3 ,ε ∈ R.
∂t
∂x
∂y

(1)

The equation and its generalizations have been used in various areas, suchas laser [4] and nonlinear optics [5] and the theory of pattern formation [2],[6].
The equation (1) admits the following traveling wave reduction [7],U 0000 + aU 00 + U 3 − U = 0, 0 :=

d
.
dZ

(2)

In our paper, we considered the equation (2) to be defined on complex
domain. By using a method of Eremenko [3], we show that all meromorphic
solutions of the traveling wave reduction (2) of (1) are elliptic or degenerate∗

Partially supported by PROCORE - France/Hong Kong joint research grant FHK29/05T and RGC grant HKU 703807P.
†
Partially supported by PROCORE - France/Hong Kong joint research grant FHK29/05T and RGC grant HKU 703807P.
‡
Partially supported by RGC grant HKU 703807P and a post-graduate studentship at
HKU.

53
1

elliptic. The method is based on the analysis using Nevanlinna theory fromcomplex analysis. We then obtain all the meromorphic solutions explicitly
by the subequation method of Conte and Musette [1]. This is an algorithmic
method, so is systematic. One of the solutions appears to be a new elliptic
solution.

References[1] R. Conte and M. Musette, Analytic solitary waves of nonintegrable equations, Phys. D 181 (2003) pp. 70–79.
[2] M.C. Cross and P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Modern Phys. 65 (1993) pp. 851–1112.
[3] A.E. Eremenko, Meromorphic traveling wave solutions of the KuramotoSivashinsky equation, J. Math. Phys., Anal. Geom. 2 (2006) pp. 278–
286.
[4] J. Lega, J.V. Moloney and A.C. Newell, Swift-Hohenberg equation for
lasers, Phys. Rev. Lett. 73 (1994) pp. 2978–2981.
[5] S. Longhi and A. Geraci, Swift-Hohenberg equation for optical parametric oscillators, Phys. Rev. A 54 (1996) pp. 4581–4584.
[6] R. Hoyle, Pattern formation: an introduction to methods (Cambridge
University Press, Cambridge, 2006).
[7] L.A. Peletier and W.C. Troy, Spatial patterns: higher order models inphysics and mechanics (Springer, 2001).

54

Free-surface flow under a sluice gateof an inclined wall from deep water
L.H. Wiryanto
Faculty of Mathematics and Natural Sciences, Bandung Institute of Technology
e-mail: leo@math.itb.ac.id
Abstract
Nonlinear solutions of free surface flow under a sluice gate are studied in thispaper. Upstream, the fluid is assumed to be infinitely in depth, and the gate makes
an angle β to the horizontal axis. Therefore, the flow emerges from the gate andproduces uniform stream far downstream. The problem is solved numerically by aboundary element method derived from an integral equation along the free surface.
An analytical function is constructed, relating to the upstream flow, so that the
integral equation is solvable. As the result, a free surface flow with smoothdetachment from the edge of the gate is obtained for relatively large upstream
Froude numbers, otherwise a free surface with back flow near the edge of the gate
is indicated, and it tends to a stagnation point for a certain Froude number.

55

Free Convection Boundary Layer Flow from a
Vertical Cone in a Porous Medium Filledwith a Nanofluid
S. Ahmad1),∗ , I. Pop2) , K. Vajravelu3)
1) School of Mathematical Sciences, Universiti Sains Malaysia, 11800 USM,
Penang, Malaysia
2) Faculty of Mathematics, University of Cluj, R-400082 Cluj, CP 253,
Romania3) Department of Mathematics, University of Central Florida, Orlando, FL
32816, USA∗ e-mail: syakilaahmad@usm.my
Abstract
A study of convective flow in porous media filled with conventional heattransfer fluids has been widely conducted in the recent years due to its wide
applications in engineering as post-accidental heat removal in nuclear reactors, solar collectors, drying processes, heat exchangers, geothermal and oil
recovery, building construction, etc. (see [1], [2], [3], [4]). It is well known
that conventional heat transfer fluids including water, oil and ethylene glycol
mixture are poor heat transfer fluids, since the thermal conductivity of these
fluids plays an important role on the heat transfer coefficient between theheat transfer medium and the heat transfer surface. In order to improve theefficiency of the fluids in the term of heat transfer, a new kind of heat transfer fluids known as nanofluid which contains a small quantity of nanosizedparticles (usually less than 100nm) that are uniformly and stably suspended
in a liquid has been introduced by Choi [5].
In this paper, a problem of steady free convection boundary layer flownear a vertical cone (with half angle σ) embedded in a saturated porous
medium filled with a nanofluid is theoretically studied using three typesof metallic and non-metallic nanoparticles such as copper (Cu), alumina(Al2 O3 ) and titania (T iO2 ) in the base fluid of water. It appears that thisproblem have been investigated by Cheng et al. [6], Yih [7] and Sohouli etal. [8] but only for the conventional heat transfer fluids. Therefore, our aim
is to extend the problem by considering a nanofluid with constant surface
heat flux. Under the assumption of using the model proposed by Tiwari
and Das [9] along with the assumptions of the Darcy-Boussinesq and boundary layer approximations, the basic partial equations along with the boundary conditions are obtained. Then, these equations are reduced to ordinary56

differential equations (ODEs) by introducing the similarity variables. The
ODEs are solved numerically using shooting method for some values of thenanoparticle volume fraction parameter ϕ. The local Nusselt number andthe temperature profiles are presented and discussed in detail. It is foundthat the local Nusselt number and the temperature profiles increase as thevolume fraction parameter increases.

References[1] D.A. Nield, A. Bejan: Convection in Porous Media, 3rd ed., Springer,
New York, 2006.
[2] D.B. Ingham, I. Pop (Eds.): Transport Phenomena in Porous Media,
Vol. III, Elsevier, Oxford, 2005.
[3] K. Vafai (Ed.): Handbook of Porous Media, 2nd ed., Taylor & Frabcis,
New York, 2005.
[4] P. Vadasz: Emerging Topics in Heat and Mass Transfer in Porous Media,
Springer, New York, 2008.
[5] S.U.S. Choi: Enhancing Thermal Conductivity of Fluids with Nanoparticles, in: D.A. Siginer, H.P. Wang (Eds.), Development and Applications of non-Newtonian Flows, ASME MD-vol. 231 and FED-vol. 66(1995), 99–105.
[6] P. Cheng, T.T. Le, I. Pop: Natural Convection of a Darcian Fluid abouta Cone, International Communications in Heat and Mass Transfer 12(1985), 705–717.
[7] K.A. Yih: Uniform Transpiration Effect on Combined Heat and Mass
Transfer by Natural Convection over a Cone in Saturated Porous Media:
Uniform Wall Temperature/Concentration or Heat/Mass Flux, International Journal of Heat and Mass Transfer 42 (1999), 3533–3537.
[8] A.R. Sohouli, M. Famouri, A. Kimiaeifar, G. Domairry: Applicationof Homotopy Analysis Method for Natural Convection of Darcian fluid
about a vertical full cone embedded in Porous Media Prescribed Surface
Heat Flux, 15 (2010), 1691–1699.
[9] R.K. Tiwari, M.K. Das: Heat Transfer Augmentation in a Two-sided
Lid-driven differentially Heated Square Cavity Utilizing Nanofluids, International Journal of Heat and Mass Transfer 50 (2007), 2002–2018.

57

Features of the Nyström Method for the
Sherman-Lauricella Equation on Piecewise
Smooth Contours
Didenko Victor1),∗ , Helsing Johan2),⋆1) Faculty of Science, University of Brunei Darussalam2) Centre for Mathematical Sciences, Lund University∗ e-mail: diviol@gmail.com⋆ e-mail: helsing@maths.lth.se
AbstractLet Γ be a simple closed contour in the complex plane C. The ShermanLauricella equation
(
)
(
)
∫
∫
τ −t
τ −t
1
1
ω(t) +
ω(τ ) d ln−
ω(τ ) d= f (t), t ∈ Γ.
2πi Γ
2πi Γ
τ −t
τ −t
(1)
plays an important role in various applications. Thus it is used in radar
imaging, in theory of viscous incompressible flows, as well as in plane elasticity and in other problems of mathematical physics [1, 2]. However, thesolution of equation (1), if it exists, is not unique and there is no analytic
formula to obtain it.
In this work we study the stability of an approximation method in thespace L2 (Γ) for equation (1) in the case where contour Γ has corner points.
More precisely, we consider the Nyström method based on composite GaussLegendre quadrature formula∫ 1u(s) ds ≈
0

n−1 ∑
d−1
∑

wp u(slp )/n,(2)

l=0 p=0

where

l + εp, l = 0, 1, . . . , n − 1; p = 0, 1, . . . , d − 1,(3)nand wp and 0 < ε0 < ε1 , . . . < εd < 1 are the weights and the Gauss-Legendre
points on the interval [0, 1].
It is worth noting that for smooth contours Γ the integral operators in theleft-hand side of (1) are compact, and the method proposed is always stable.
On the other hand, if Γ possesses angular points τk , k = 1, 2, . . . , m, the stability of the method under consideration is connected to the invertibility ofcertain operators Aτk from an algebra of Toeplitz operators. The operators
slp =

58

Aτk depend on the above parameters εp and on the magnitude of the corresponding angles at the points τk . These operators Aτk arise if one applies
quadrature formulas similar to formulas (2)-(3) to Mellin operators with thekernels defined by the kernels of the integral operators in equation (1). They
have a complicated structure and, at the moment, there is no analytic tools
to study their invertibility. However, each such operator can be associatedwith the corresponding Nyström method for the Sherman-Laurichella equation on a special model curve Γ0 with only one corner point. On the other
hand, the stability of the method in those model situations is connected tothe behaviour of condition numbers of systems of algebraic equations andit can be tested numerically. Thus, considering the sequences of the condition numbers we obtain certain information about the invertibility of theoperators Aτk .

References[1] S. G. Mikhlin: Integral Equations and Their Applications to Certain
Problems in Mechanics, Mathematical Physics and Technology, 2nd ed.,Pergamon Press, 1964.
[2] N. I. Muskhelishvili: Fundamental Problems in the Theory of Elasticity,
Nauka, 1966.

59

A New Evolutionary Algorithm for Groundto Flight Data Correlation
Changtong Luo1),∗ , Shao-Liang Zhang2) , Zonglin Jiang1)
1) Institute of Mechanics, Chinese Academy of Sciences
2) Department of Computational Science and Engineering, Nagoya
University∗ e-mail: luo@imech.ac.cn
Abstract
In hypersonic vehicle design, it is very difficult/expensive to simulate
the environment in flight. Therefore, it is of key important to predict the
aerodynamic coefficients of the flight vehicle using wind-tunnel and flighttest data. The task is called ground to flight data correlation, also shorten
as ground/flight correlation. A number of methods have been applied to
ground/flight correlation, including extrapolation [1], least squares regression [2], artificial neural network [3, 4] and maximum likelihood method [5].
However, all these methods are based on a certain model structure and optimize the parameters in the model. This is, in many cases, unreliable. To
find an optimal model structure will definitely improve the predicted results.
In this study, we suggest searching the best model structure from a functional space F with collected data {{xi1 , xi2 , · · · , xin ; yi }}N
i=1 , which can be
described as f ∗ = arg min
f ∈F

N
X

kf (xi1 , xi2 , · · · , xin ) − yi k .

iGenetic programming [6] is a candidate for this optimization problem.
However, its tree-based representation makes it difficult (although not impossible) to be implemented in general-purpose programming languages suchas C/C++ and Fortran. To overcome this difficulty, M. O’Neill and C. Ryan
proposed a grammatical evolution (GE) [7] in 2001. However, GE is still notso easy to use. First, the implementation of GE is complicated because it
needs an additional function parser for the encoding and decoding process.
Next, the incomplete mapping and extra codons problems [7, 8] are common
but difficult to handle.
In this paper, a new evolutionary algorithm, parse-matrix evolution (PME),
for flight data correlation is proposed. A chromosome in PME is a parsematrix with integer entries. The mapping process from a chromosome to its
model function is based on a mapping table. PME can easily be implemented
in any programming language and free to control. Furthermore, it does notneed any additional function parsing process. Numerical results show that
PME can solve the symbolic regression problems effectively (see Table 1).
60

Table 1: Test models and performance of PME
no. Dim Target model1
1x2 − sin x
2
1
sin x + 2x
31
cos x2 − x
41
sin x + cos x − x
52x1 + 2x2
62(x1 + x2 )/x272sin(x1 + x2 )2cos(x1 + x2 ) − x282ln(x1 + x2 )9102
x1 − ex1 +x2
2
ln(x1 + x2 ) + sin(x1 + x2 )11122
sin(x21 − x2 ) + ln(x1 + x2 )3
x1 + x2 − x3
133
x1 x2 − x3
14153
sin(x1 x2 ) + x3

Domain[1, 3][1, 3][1, 3][1, 3][1, 3]2[1, 3]2[1, 3]2[1, 3]2[1, 3]2[1, 3]2[1, 3]2[1, 3]2[1, 3]3[1, 3]3[1, 3]3

no.samples ave no.eval ave no.start
8
3230
4.8
8
2071
1.1
8
2603
2.3
8
3646
7.4
16
1047
1.0
16
1701
1.3
16
1537
1.0
16
2798
1.5
16
1584
1.2
16
3593
1.8
16
4911
5.7
16
7637
13.2
64
4505
1.9
64
3164
2.7
64
3302
11.6

References[1] J.B. Peterson, Jr., M.J. Mann, R.B. Sorrells III, W.C. Sawyer, and D.E.
Fuller: Extrapolation of Wind-Tunnel Data to Full-Scale Conditions, NASA
TP-1515, 1980.
[2] E.A. Morelli, R. DeLoach: Wind Tunnel Database Development Using Modern Experiment Design and Multivariate Orthogonal Functions, AIAA paper,
AIAA 2003-0653, 2003.
[3] M. Norgaard, C. Jorgensen, J. Ross: Neural Network Prediction of New
Aircraft Design Coefficients, NASA Technical Memorandum 112197, 1997.
[4] T. Rajkumar, Jorge Bardina: Prediction of Aerodynamic Coefficients using
Neural Networks for Sparse Data, in: Proc. 15th Int’l Florida Artificial Intelligence Research Society Conf., 242-246, 2002.
[5] J.-H. Lee, E.T. Kim, B.-H. Chang, I.-H. Hwang, D.-S. Lee, The Accuracy of
the Flight Derivative Estimates Derived from Flight Data, World Academyof Science, Eng. Technol., 58 (2009), 843-849.
[6] J. R. Koza, Genetic programming: on the programming of computers bymeans of natural selection. Cambridge, MA: MIT Press, 1992.
[7] M. O’Neill and C. Ryan, Grammatical evolution, IEEE Trans. Evol. Comput.
5 (2001), 349-358.
[8] M. O’Neill, A. Brabazon, Grammatical swarm: the generation of programsby social programming. Nat. Comput. 5 (2006), 443-462.

61

Ocean wave scattering within or over a
porous medium
Swaroop Nandan Bora1),∗ , Santu Das2)
1) Department of Mathematics, Indian Institute of Tech. Guwahati, India
2) Department of Mathematics, Indian Institute of Tech. Guwahati, India∗ e-mail: swaroop@iitg.ernet.in
Abstract
Water wave scattering within a porous medium on a sea-bed is consideredbased on linear wave theory. Using the model of wave-induced flow within aporous medium and Galerkin eigenfunction expansions, refraction-diffraction
equations for surface waves are derived whenever needed. The reflection andtransmission of waves are considered mainly for two cases: (i) the porous
medium is obstructed on the right by a solid wall which does not allow any
transmission for the case of horizontal bottom, (ii) and when the bottom
is approximated by a line of constant slope or by a raised horizontal bottom. Further, by introducing a non-porous medium again after the porous
medium, subsequent scattering is investigated. The fluid region is divided
into two or three sub-regions according to the points of discontinuity. For the
first case, the wave motion inside the porous structure tends to decay as itpropagates through the pores. The wave, on encountering the solid wall, getsreflected into the porous structure. For the second case, refraction-diffraction
equations for surface waves are derived for the porous region. The scattering
is now due to both the porous structure and the structure at the bottom.
A boundary value is set up for a flat sea-bed and the analytical solution obtained and Galerkin eigenfunction method is introduced in order toexpand the spatially dependent potential in terms of depth-dependent functions. The oblique incident wave propagates at an angle θ to the x-direction.
By appropriate matching of the properties such as velocity potential, velocity
components of two nearby regions, the reflection phenomenon is studied.
For a water depth of h, which may be uniform (even bed) or non-uniform(uneven bed), the potential φ in each region satisfies Laplace’s equation∇2 φ = 0, −h ≤ z ≤ 0.

(1)

The free surface boundary condition will be∂φ− µφ = 0, at z = 0, in the non-porous medium, and,
∂z

62

(2)

∂φ
− iµR = 0, at z = 0, in the porous medium.
(3)∂zwhere µ = σ 2 /g with σ as the wave frequency, g the gravity constant, R =f − iS the impedance of the porous medium with f as friction factor and S
as inertial coefficient. Considering the depth-dependent part of the potential
cosh kn (hj + z)
as Zj,n (z) =, n = 0, 1, 2, . . . with j = 1, 2, . . . denoting thecosh kn hj
region number and the wave number kn satisfying the complex dispersion
relation iRµ = kn tanh kn hj .
In the first region prior to the location of the porous structure, the velocity
potential isφ(x, z) = (eiK1,0 x + R0 e−iK1,0 x )Z1,0 (h1 , z) +∞X

Rn e−iK1,n x Z1,n (h1 , z), (4)

n=1

q

2
where K1,n = k1,n− λ2 , with λ = k0 sin θ and R0 the reflection coefficient.
The potentials in the other regions can be found by looking at the relevant
physics. Ultimately the reflection coefficient can be obtained with the helpof coefficients through the orthogonal functions Zn (z). Some values for the
relevant parameters are assumed and the reflection is studied for various
depths, porosity factor etc.
When the wave is allowed to propagate over an uneven porous sea-bed,
perturbation technique is used to reduce the original problem to a simpler one
for the first order correction of the potentials. The solution of this problem
is obtained by an approximate use of Green’s integral theorem. Sinusoidal
ripples is considered as an example of the sea-bed. The reflection and transmission coefficients are evaluated against various values of parameters and it
is observed that very few ripples are needed to generate significant reflection.

References[1] Martha SC, Bora SN, Chakrabarti A. Oblique water wave scattering bysmall undulations on a porous sea-bed. Applied Ocean Research, 2007;
29: 86-90.
[2] Dalrymple RA, Losada MA, Martin PA. Reflection and transmission
from porous structures under oblique wave attack. J. Fluid Mech., 1991;
224: 625-644.
[3] Zhu S. Water waves within a porous medium on an undulating bed.
Coastal Engineering, 2001; 42: 87-101.

63

Spatially localized traveling structures and the asymptotic
behavior in binary fluid convection
1),*

Yasumasa Nishiura2)2), Takeshi Watanabe , Makoto Iima

1) Research Institute for Electronic Science, Hokkaido University2) Department of Mathematical and Life Sciences, Hiroshima University* e-mail: nishiura@nsc.es.hokudai.ac.jp
Abstract
We study spontaneous pattern formation and its asymptotic behavior in binary
fluid flow driven by a temperature gradient [1]. When the conductive state is
unstable and the size of the domain is large enough, finitely many number ofspatially localized time-periodic traveling structures (PTP) [3] are generated
spontaneously in the conductive state. They collide and merge via strong collisions,and finally, they are arranged at non-uniform intervals and move in the same
direction as in Fig.2. Although there are many other stable patterns, it is remarkablethat PTPs end up the winner for survival via interactions as in Fig.1. We found thatthe role of PTP solutions and their strong interactions (collision) [2] are importantto characterize the asymptotic state and hence, pulse-pulse interactions areinvestigated in detail.

Fig.1: Collision between stable stationary convective pattern (middle) and PTP(from right). They merge into one body, but PTP persists after strong interaction.

64

Fig.2: Time evolution for two generic initial data. A point disturbance is added at t
= 0, x = 500 (right), and small random disturbances are added in the whole domain(left). Time interval of the visualization is t = 250.

References[1] P. Kolodner, Collisions between pulses of traveling-wave convection, Phys.
Rev. A 44, (1991) 6466–6479.
[2] Y. Nishiura, T. Teramoto, and K. Ueda, Scattering of traveling spots indissipative system, Chaos 15, (2005) 047509.
[3] T. Watanabe, K. Toyabe, M. Iima and Y. Nishiura Time-periodic travelingsolutions of localized convection cells in binary fluid mixture, Theoreticaland Applied Mechanics Japan, 59, (2010), 211–219.
65

1

Adaptive Mixed Finite Element Methods andNew Discontinuous Galerkin methods
Eun-Jae Park1),∗ ,1) Department of Computational Science and Engineering,
Yonsei University, Seoul 120-749, Korea∗ e-mail: ejpark@yonsei.ac.kr
Abstract
This extended abstract deals with locally conservative methods and consists of two parts: ﬁrst part on adaptive mixed ﬁnite element methods andsecond on new discontinuous Galerkin methods.
First part is based on joint work with Dongho Kim [2]. We study adaptive mixed ﬁnite element discretizations for nonlinear elliptic problems. The
mixed ﬁnite element method has two important features; it conserves themass locally and produces accurate ﬂux even for highly nonhomogeneous
media with large jumps in the physical properties. In this paper, we take anapproach based on the Brezzi-Rappaz-Raviart framework rather than using
the Brouwer ﬁxed point theory. Compared with existing literature several
improvements are made. First, the Raviart-Thomas ﬁnite element space of
all orders including the lowest order case is treated. Existence and uniqueness of the approximate solution to the model problem is proved. Next, wepoint out that compared with the work [5], the required H 5/2 (Ω) regularityof the solution weakens to H 2+ϵ (Ω).
Next, we drive an optimal order a priori error estimate measured in
the Lm (Ω)-norm in the framework of Brezzi-Rappaz-Raviart. The BrezziRappaz-Raviart theory requires of the C 1,1 -regularity of the nonlinear function; we extend the theory to allow more general C 1,α class of the function
for 0 < α ≤ 1.
Finally, we derive reliable and eﬃcient a posteriori error estimators inm
L (Ω) for the error control of our approximation to the nonlinear problem
under consideration. To the best of the authors’ knowledge, this is the ﬁrst
result obtaining a posteriori error estimator measured in the Lm (Ω)-norm
which is reliable and eﬃcient for a mixed approximation of nonlinear elliptic
problems.
Second part is based on joint work with Youngmok Jeon [4]. A newclass of hybrid discontinuous Galerkin methods is introduced and analyzed
for second-order elliptic equations. Recently, the hybridizable discontinuous66

2
Galerkin (HDG) method is introduced and developed by Cockburn and his
collaborators [1]. The main feature of the HDG method is that their approximate solutions can be expressed in an element-by-element fashion in terms ofan approximate trace ubh satisfying a global weak formulation, which reduces
globally coupled degrees of freedom dramatically. In this sense, our method
can be understood as a HDG method.
On the other hand, our method can be viewed as a generalization of socalled cell boundary element (CBE) by Jeon and his coauthors [3], therefore,
the method is named as GCBE. The CBE method is a ﬂux conserving method
without a dual partition, but it has the same disadvantages as the ﬁnite
volume method, namely that high-order accurate methods are not easy toconstruct [3]. The motivation of the GCBE method is to construct a ﬂux
conserving method with high order approximation property. More details
can be found in the recent paper [4].

References[1] B. Cockburn, J. Gopalakrishnan and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkinmethods for second order elliptic problems, SIAM J. Numer. Anal. 47(2009), pp. 1319–1365.
[2] D. Kim and E.-J. Park, A Priori and A Posteriori Analysis of Mixed Finite Element Methods for Nonlinear Elliptic Equations. SIAM J. Numer.
Anal. Volume 48 (2010), no. 3, pp. 1186-1207.
[3] Y. Jeon, E.-J. Park, Nonconforming cell boundary element methods for
elliptic problems on triangular mesh, Appl. Numer. Math. 58 (6) (2008)
800–814.
[4] Y. Jeon and E.-J. Park, A hybrid discontinuous Galerkin method for
elliptic problems, SIAM J. Numer. Anal. Volume 48 (2010), no. 5, pp.
1968-1983[5] F. A. Milner, E.-J. Park, Mixed finite element methods for HamiltonJacobi-Bellman-type equations, IMA J. Numer. Anal. Vol.16 (1996), pp.
399–412.

67

Numerical Investigation of Mixed Convection
Flow over a Horizontal Circular Cylinderand a Solid Sphere in Nanofluids
Roslinda Nazar
School of Mathematical Sciences, Faculty of Science & Technology
Universiti Kebangsaan Malaysia
e-mail: rmn@ukm.my
Abstract
The recent discovery of nanofluids by Choi [1], which is a new kind of fluid
suspension consisting of uniformly dispersed and suspended nanometer-sized(10-50 nm) particles and fibers in the base fluid, marks the next approachas a cooling technology. This is because nanofluids exerting nano effects
have special features. One such characteristic of nanofluids is the anomalous
high thermal conductivity at very low concentration of nanoparticles and theconsiderable enhancement of forced convective heat transfer, as it was pointed
out by many researchers. Nanofluids usually contain the nanoparticles suchas metals, oxides, or carbon nanotubes, whereby these nanoparticles haveunique chemical and physical properties.
In this study, the steady mixed convection boundary layer flow of nanofluids, formed by the dilution of nanoparticles, past a horizontal circular cylinderand about a solid sphere with constant surface temperature is investigatednumerically for both cases of assisting and opposing flows. The nanofluid
equations model proposed by Tiwari and Das [2] is used for the analysis.
The resulting system of nonlinear partial differential equations is solved using an implicit finite-difference scheme known as the Keller-box method. The
solutions for the flow and heat transfer characteristics are evaluated numerically for various values of the parameters, namely the nanoparticle volume
fraction ϕ and the mixed convection parameter λ at Prandtl number Pr =0.7, 1, 6.2 and 7. Three different types of nanoparticles, namely Al2 O3 , Cu
and T iO2 (by using water-based fluid with Pr = 6.2) have been consideredin this study and representative results for the local skin friction coefficientand the local heat transfer coefficient or local Nusselt number have been obtained for the following range of nanoparticle volume fraction: 0 ≤ ϕ ≤ 0.2(see Abu-Nada and Oztop [3]). We have used data related to thermophysicalproperties of fluid and nanoparticles as presented in [3] to compute each caseof nanofluid.
In order to verify the accuracy of the present method, the values of thelocal skin friction coefficient and the local Nusselt number are compared68with those reported by Merkin [4] for ϕ = 0 (regular Newtonian fluid) at
Pr = 1 and various values of x and λ, for the horizontal circular cylindercase; while the values of the local skin friction coefficient and the local heattransfer coefficient are compared with those reported by Nazar et al. [5]for ϕ = 0 at Pr = 0.7 and Pr = 7 and various values of x and λ, for thesolid sphere case. The comparison of the previous results obtained in [4, 5]
and the present ones are found to be in excellent agreement. It is foundthat for each particular nanoparticle, as the nanoparticle volume fraction
ϕ increases, the skin friction coefficient and the heat transfer rate at thesurface also increase, and it also leads to the increase of the value of themixed convection parameter λ which first gives no separation and the valueof λ below which a boundary layer solution does not exist. It is also foundthat an increase in the Prandtl number Pr leads to an increase of the local
skin friction coefficient and the local Nusselt number. Finally, it is found thatthe velocity and temperature profiles satisfy the far field boundary conditions
asymptotically, and as such, this support the validity of the numerical results
obtained. The results obtained can be used to explain the characteristics andapplications of nanofluids, which are widely used as coolants, lubricants, heatexchangers and micro-channel heat sinks.

References[1] S. Choi: Enhancing thermal conductivity of fluids with nanoparticles,
ASME International Mechanical Engineering Congress and Exposition,
FED 231/MD, 66 (1995) 99–105.
[2] R.K. Tiwari, M.K Das: Heat transfer augmentation in a two-sided liddriven differentially heated square cavity utilizing nanofluids, International Journal of Heat and Mass Transfer, 50 (2007), 2002–2018.
[3] E. Abu-Nada, H.F. Oztop: Effects of inclination angle on natural convection in enclosures filled with Cu-water nanofluids, International Journalof Heat and Mass Transfer, 30 (2009), 669–678.
[4] J.H. Merkin: Mixed convection from a horizontal circular cylinder, International Journal of Heat and Mass Transfer, 20 (1977), 73–77.
[5] R. Nazar, N. Amin, I. Pop: Mixed convection boundary layer flow about
an isothermal sphere in a micropolar fluid, International Journal of Thermal Sciences, 42 (2002) 283–293.

69

Lamellar pattern and two-phase shallow
water equations
Hidenori Yasuda1)∗
1) Department of Mathematics, Josai University∗ e-mail: yasuda@math.josai.ac.jp
Abstract
We propose a numerical method to simulate the phase separation of block
copolymer in one-dimensional thin liquid ﬁlm.
Phase separation of polymer blend: a mixture of homopolymers, ﬁnallybecomes a single circle [1]. However, phase separation of block copolymer:
chemically connected two homopolymers, is qualitatively diﬀerent. In onedimensional system, periodic microscopic scale pattern called lamellar pattern appears.
Phenomena of phase separation of copolymer is modeled by two-phaseshallow water equations and Ohta-Kawasaki potential. Ohta-Kawasaki potential governs the phase separation of block copolymer [2]. However, somediﬃculty appears in the computation of Ohta-Kawasaki potential. A modiﬁcation of the potential is proposed to circumvent the diﬃculty. Twophase shallow water equations are discretized by an invariant ﬁnite diﬀerence
method. And the potential is integrated by the trapezoidal rule. Lamellar
pattern and wavy surface developed using the proposed method [3].
Models and numerical method
1. Two-phase shallow water equations
Two-phase ﬂow in the liquid ﬁlm is described by two-phase shallow waterequations,∂ (αd h)∂+(αd hud ) = 0,
(1)
∂t
∂x
∂
∂ (αc h)+(αc huc ) = 0,
(2)
∂t
∂x
∂ud
∂
1 ∂
∂2
+ ud ud = −(ρm gh) + ν 2 ud + fd ,
(3)
∂t
∂x
ρd ∂x
∂x
∂
1 ∂
∂2
∂uc
+ uc u c = −(ρm gh) + ν 2 uc + fc ,
∂t
∂x
ρc ∂x
∂x

(4)

αd + αc = 1.

(5)

70

Here, α is the volume fraction, u the velocity, h surface height, ρ density,d is the suﬃx of minor phase, c is the suﬃx of major phase, and ρm is thedensity of mixture:
Two-phase shallow water equations is discretized by an invariant ﬁnite diﬀerence scheme which highly resolve the interface region of phase separations [1].
2. Ohta-Kawasaki potential
Phase separation of block copolymer is governed by Ohta-Kawasaki potential [2],F {ψ} = Fs {ψ} + Fl {ψ}(6)ψ = αd − αc
Fs {ψ} =Fl {ψ} =

∫

[

a
b
γ
dx − ψ 2 + ψ 4 + (∇ψ)2
2
4
2

]

∫
(
)(
)
ε∫
dx dx0 G (x, x0 ) ψ (x) − ψ̄ ψ (x0 ) − ψ̄2−∆G (x, x0 ) = δ (x − x0 ) .

(7)
(8)
(9)
(10)

Here, F means free energy of block copolymer; ψ is the order parameter: ψ̄
is the average of ψ. G is the fundamental solution of Laplace equation.
Ohta and Kawasaki potential has some diﬃculty to compute numerically.
For lamellar patterns, we propose a modiﬁcation to circumvent the diﬃculty.
Numerical integration of the integration of the potential is performed by thetrapezoidal rule.
In simulations using the proposed method, lamellar pattern and wavy surface
consistent with the speciﬁed period developed [3].

References[1] H. Yasuda, Two-phase shallow water equations and phase separation inthin immiscible liquid ﬁlms, J. Sci. Comput. 43 (2010) 471-487.
[2] T.Ohta, K. Kawasaki, Equilibrium morphology of block copolymer melts,
Macromolecules, 19 (1986) 2621-2632.
[3] H. Yasuda, Simulation of micro-phase separation of block copolymer in
one-dimensional thin liquid ﬁlms, submitted.

71

IMPLEMENTING THE FINITE VOLUME METHOD ON A RECTANGULAR GRID

Tanay Milind Deshpande*
Undergraduate, B.E.(Hons.) Mechanical Engineering
Birla Institute of Technology and Science (BITS) – Pilani K. K. Birla Goa Campus
Goa, India*tanay.d255@gmail.com

Abstract
The Finite Volume Method (FVM) is commercially the most widely used CFD technique toconvert the non-linear partial differential equations that govern fluid flow into a linear
algebraic system of equations. This paper formulates an algorithm to solve steady-state,laminar and incompressible fluid flow equations with a structured collocated gridarrangement in a rectangular domain and describes its implementation using C++. The chief
aim of this algorithm is to considerably simplify and optimize existing FVM techniques suchas the Semi Implicit Method for Pressure Linked Equations (SIMPLE) algorithm (Patankar,1972) for low Peclet number problems in regular geometries. After setting boundary
conditions for driven-cavity flow as an instance, the Central Differencing scheme and thedivergence theorem are invoked. The Navier-Stokes equations are discretized at each cell in
the grid into a system of coefficients of neighbouring cell velocities. The algebraic equations
representing each line in the grid are then fed into a solver subroutine that implements the
Tridiagonal Matrix Algorithm. After two Alternating Direction Implicit traversals for velocities,
the continuity equation is discretized in a similar manner for evaluating pressure correction
terms which update the pressure and velocities at every node for the next iteration. At thiscrucial juncture where the SIMPLE takes contributions of pressure errors at neighbouring
cells to be zero while calculating that at a single cell, this algorithm departs from theprocedure and instead accommodates for neighbouring pressure errors.
References1. Versteeg H., Malalasekera W., ‘Introduction to Computational Fluid Dynamics: The
Finite Volume Method’, 1995, Longman Scientific & Technical
2. Anderson J., ‘Computational Fluid Dynamics: The Basics with Applications’, 1995,
McGraw-Hill
3. Patankar S., ‘Numerical Fluid Flow and Heat Transfer’, Hemisphere Publishing
Corporation, 1980
4. Ferziger J., Peric M., ‘Computational Methods for Fluid Dynamics, Springer-Verlag,
2001
5. The SIMPLE Algorithm, Patakar S., Spalding B., 1972 & The Tridiagonal Matrix
Algorithm, Thomas, 1949

72

Accurate and Rigorous Exponential
Algorithm in Round to Nearest
Naoya Yamanaka1),∗ , Masahide Kashiwagi2) , Shin’ichi Oishi2)
1) Research Institute for Science and Engineering, Waseda University
2) Faculty of Science and Engineering, Waseda University∗ e-mail: naoya yamanaka@suou.waseda.jp
Abstract
This talk is concerned with numerical algorithms retaining high reliability,
high accuracy and high portability. In this talk, an algorithm with highreliability means a numerical algorithm which outputs a mathematicallyrigorous result. An algorithm with high accuracy represents an algorithm
which returns a result with high accuracy. Furthermore, an algorithm withhigh portability indicates an algorithm which calculates a result without
relying on any numerical environment. In this talk, numerical exponential
algorithm retaining high reliability, high accuracy and high portability is
discussed.
Todays libraries for the approximation of elementary functions are very
fast and the results are mostly of very high accuracy [1, 2]. The achieved
accuracy does not exceed one or two ulp for almost all input arguments; however, there is no proof for that. In this talk, we discuss a reliable, accurate,
portable implementation of the exponential function. We intend to utilize
the marvelous accuracy by a table approach. Proposed algorithm delivers
rigorous bounds for the result for all floating point input arguments. The
order of evaluation of the formula is carefully chosen to diminish accumulation of rounding errors. As a result we obtain the relative accuracy of the
exponential function value is better than 1ulp.
Throughout this talk, we assume floating-point arithmetic adhering toIEEE standard 754-1985. IEEE standard 754-1985 is one of a technicalstandard established by the IEEE and the most widely-used standard for
floating point computation, followed by many hardware and software implementations. Proposed algorithm is based on computer interval arithmetic
adhering to the standard. Interval arithmetic is a numerical method developed by mathematicians since the 1950s and 1960s as an approach to puttingbounds on rounding errors and measurement errors in mathematical computation and thus developing numerical methods that yield reliable results.
It represents each value as a range of possibilities. The interval arithmetic
on floating-point arithmetic with changing rounding mode is widely used.
73

Rounding mode is one of defined parameter in IEEE standard 754-1985. The
standard defines four rounding algorithms: rounding to nearest, round toward zero, round toward plus-infinity, round toward minus-infinity. However,changing rounding mode takes some computational costs, and the commands
to change the mode vary widely depending on the numerical environment.
Besides, numerical environments which do not have the commands to change
the mode exist. For these problems, proposed algorithm only in rounding tonearest is proposed.

References[1] J. M. Muller: Elementary Functions, Algorithm and Implementation,
Second Edition, Birkhäuser, 2006.
[2] S. M. Rump: Rigorous and portable standard functions. BIT Numerical
Mathematics, 41(3):540–562, 2001.

74

An integrated manufacturing system fortime-varying demand process
M.OmarInstitute of Mathematical Sciences, University of Malaya, Kuala Lumpur
e-mail: mohd@um.edu.my
Abstract
This paper considers a manufacturing system in which a single-manufacturer
procures raw material from a single-supplier at multiple installments andprocess them to make finished products and deliver to a single-buyer at multiple shipments to satisfy a time-varying demand rate. In this system, themanufacturer must deliver the products in small quantities to minimize thebuyer’s holding cost, and accept the supply of small quantities of raw material to minimize its own holding cost. We develop a mathematical modelfor this problem and illustrate the effectiveness of the model with numericalexamples.

75Numerical analysis of a stationary transport
equation by upwind scheme
Takanori Katsura
Graduate School of Informatics, Kyoto University
e-mail: katsura@acs.i.kyoto-u.ac.jp
Abstract
The aim of this study is to simulate the stationary transport equation
numerically as a fundamental study of the diﬀuse optical tomography(DOT).
DOT is a new biomedical imaging technique, which is hoped to come after
X-ray CT and MRI, and it has been extensively studied for practical use
in medical sciences. It is considered that light propagation in biomedical
tissue is described as photon density by the transport equation, and DOT
is mathematically formulated as an inverse problem of the equation [1, 2].
This mathematical model of DOT is commonly accepted but should be more
veriﬁed in detail to put DOT to practical use. To this end, we propose a
new numerical scheme based on the ﬁnite diﬀerence method which is helpful
for a fundamental study of the equation. The comparison of the numericalsimulations and the experimental measurements will give us some suggestions
for the mathematical modeling.
In this presentation, we consider the following boundary value problem
of stationary transport equation:∫ξ ·∇x u(x, ξ) + µt (x)u(x, ξ) = µs (x)p(ξ · ξ 0 )u(x, ξ 0 ) dσξ0 + f (x) (x, ξ) ∈ Q
S10u(x, ξ) = α u(x, ξ ) + q(x, ξ)(x, ξ) ∈ Γ− .

Here Ω is a square domain (0, 1) × (0, 1) and n(x) is the outward unit normal
vector on ∂Ω, Γ+ is the outward boundary, Γ− is the inward boundary;
Γ± := {(x, ξ) ∈ ∂Ω×S 1 |ξ·n(x) < 0}, α is reﬂectance, ξ 0 = ξ−2 (ξ·n(x)) n(x).
To discretize the transport equation, we apply the trapezoidal formula
to the integral term, and the upwind diﬀerence to ∇x . Let K be the division number of velocity and h be the the spatial mesh size. We set ξk :=, sin 2π
) ∈ S 1 , xi,j := (ih, jh) ∈ Ω, s(k, i) := sgn(ξk,i ). Then we pro(cos 2π
KK
pose the following upwind scheme (1)−(4):ξk,1ui,j,k − ui−s(k,1),j,kui,j,k − ui,j−s(k,2),k+ ξk,2+ µt (xi,j )ui,j,ks(k, 1)hs(k, 2)h
(1)
K−1
2π ∑
= µs (xi,j )p(ξk · ξl )ui,j,l + f (xi,j ) (xi,j , ξk ) ∈ Q
K l=0
76

ξk,1

ui±1,j,k − ui±1,j−s(k,2),kui,j,k − ui±1,j,k+ ξk,2+ µt (xi±1,j )ui±1,j,k∓hs(k, 2)h
K−1
2π ∑
(2)
= µs (xi±1,j )p(ξk · ξl )ui±1,j,l + f (xi±1,j )K l=0(xi,j , ξk ) ∈ Γ+ (x1 = 0, x1 = 1)
ui,j±1,k − ui−s(k,1),j±1,kui,j,k − ui,j±1,k
ξk,1+ ξk,2+ µt (xi,j±1 )ui,j±1,ks(k, 1)h
∓h
K−1
2π ∑
(3)
= µs (xi,j±1 )p(ξk · ξl )ui,j±1,l + f (xi,j±1 )K l=0(xi,j , ξk ) ∈ Γ+ (x2 = 0, x2 = 1)
ui,j,k = α ui,j,l + q(xi,j , ξk )(4)
ξl = ξk − 2 (ξk · n(xi,j )) n(xi,j ), (xi,j , ξk ) ∈ Γ− .
We prove the convergency of this upwind scheme. Moreover we propose aiterative scheme based on the upwind scheme and prove the convergency.
The following ﬁgures are numerical examples of the light intensity U (x) by
our proposed scheme. Light scattering is described by the scattering kernelp(ξ, ξ 0 ) and the anisotropy factor g controls how strongly the scattering is
forward-peaked; g = 0, isotropic scattering; g = 0.9, strongly forward-peaked
scattering.

g = 0.0

g = 0.5

g = 0.9

Fig. 1 Numerical examples of the light intensity U (x)References[1] S. R. Arridge, M. Schwiger: Image reconstruction in optical tomography
Philosophical Transactions of the Royal Society B, 352 (1997), 717-726[2] A. D. Klose, E. W. Larsen: Light transport in biological tissue based on
the simpliﬁed spherical harmonic equations, Journal of Computational
Physics, 220 (2006), 441-470

77

Stable Traffic Jam Occurs from Nonlinear
Saturation
Akiyasu Tomoeda1,2),∗ , Ryosuke Nishi3,4) ,Kazumichi Ohtsuka5,6) , Katsuhiro Nishinari3,5)1) Meiji Institute for Advanced Study of Mathematical Sciences, Meiji
University
2) JST, CREST
3) School of Engineering, The University of Tokyo
4) Japan Society for the Promotion of Science5) Research Center for Advanced Science and Technology, The Universityof Tokyo6) Economics and Social Research Institute, Cabinet office, Government ofJapan∗ e-mail: atom@isc.meiji.ac.jp
Abstract
Among various kinds of jamming phenomena, a traffic jam of vehicles isa very familiar phenomenon and causes several losses in our daily life such as
decreasing efficiency of transportation, waste of energy, serious environmental
degradation, etc. In particular, the dynamics of traffic flow on highway has
attracted many researchers and has been investigated as a non-equilibrium
system of interacting particles for the last few decades [1, 2] . The dynamics
of traffic flow is often treated as an effectively one-dimensional compressible
fluid by focusing on the collective behavior of vehicles, so-called macroscopic
models [3, 4] , which consist of two equations: conservation law and motion equation of traffic vehicles. In previous fluid models, one does not haveany choice to introduce the diffusion term into the models such as Kernerand Konhäuser model [4] , in order to represent the stabilized density wave,
which indicates the formation of traffic jam. However, it emerges as a serious problem that some vehicles move backward even under heavy traffic,
as mentioned in [5] . Thus, unfortunately we would have to conclude that
traffic models which include the diffusion term are not reasonable for the realistic expression of traffic flow. Given these factors, we suppose that traffic
jam forms as a result of the plateaued growth of small perturbation by thenonlinear saturation effect.
In this contribution, we have proposed a new compressible fluid model forthe one-dimensional traffic flow taking into account a variation of the reaction
time of drivers, which is based on the actual measurements. Our new model78

is a generalization of the Payne model [3] by introducing a density-dependentfunction of reaction time as a reasonable assumption, and given by∂ρ∂+(ρv) = 0,∂t ∂x
1 dVopt (ρ) ∂ρ
∂v
∂v
1 
Vopt (ρ) − v +
+v
=
,
∂t
∂x
τ (ρ)
2ρτ (ρ) dρ ∂x

(1)
(2)

where ρ(x, t) and v(x, t) correspond to the spatial vehicle density and theaverage velocity at position x and time t, respectively. τ (ρ) and Vopt (ρ)
are the reaction time of drivers and the optimal velocity function under thecurrent density ρ, respectively.
The linear stability analysis of this new model shows the instability ofa free flow state caused by the enhancement of fluctuations around a certain critical density of vehicles, which is observed in real traffic flow [6] .
Moreover, the condition of the nonlinear saturation of density against small
perturbation is theoretically derived from the reductive perturbation method.

References[1] D. Chowdhury, L. Santen and A. Schadschneider: Statistical physics
of vehicular traffic and some related systems, Phys. Rep. 329 (2000),
199–329.
[2] D. Helbing: Traffic and related self-driven many-particle systems, Rev.
Mod. Phys. 73 (2001), 1067–1141.
[3] H. J. Payne: Mathematical Models of Public Systems 1, ed G. A. Bekey(La Jolla,CA: Simulation Council) (1971), 51–61.
[4] B. S. Kerner and P. Konhäuser: Cluster effect in initially homogeneous
traffic flow, Phys. Rev. E 48 (1993), R2355–R2338.
[5] C. F. Daganzo: Requiem for second-order fluid approximation of traffic
flow, Trans. Res. B 29B (1995), 277-286.
[6] Yuki Sugiyama, et al.: Traffic jams without bottlenecks – experimental
evidence for the physical mechanism of the formation of a jam, New J.
Phys. 10 (2008), 033001.

79

Lagrangian Skeleton on the Gulf of Mexico
Oil Spill
C. Mendoza1),∗ , A.M. Mancho2)
1) Universidad Politécnica de Madrid. ETSI Navales.
Av. Arco de la Victoria s/n. 28040 Madrid, Spain
2) Instituto de Ciencias Matemáticas. CSIC-UAM-UC3M-UCM.
Serrano 121. 28006 Madrid, Spain∗ e-mail: carolina.mendoza@upm.es
Abstract
The understanding of the circulation of ocean currents, the exchange of
CO2 between atmosphere and oceans, and the influence of the oceans on the
distribution of heat on a global scale is key to our ability to predict and assess
the future evolution of climate [1, 2]. Global climate change is affecting sea
breathing through mechanisms not yet understood [3]. The ocean is important in the regulation of heat and moisture fluxes, and oceanic physical andbio-geochemical processes are major regulators of natural greenhouse gases.
Understanding how oceans mix their waters is key to provide sound forecasts on the climate [1]. Global change also affects marine biodiversity andthreatens the survival of ecosystems and exploitable resources. To predict
not only the effects of global change on the oceans, but also the responsetime of climate feedback requires to improve detection systems and to open
new lines of research.
We use a novel Lagrangian descriptor (function M, introduced in [4, 5]).
It is based on the measure of the arclength of particle trajectories on theocean surface at a given time. In [6, 7, 8] this technique has been proven tobe successful for characterizing the Kuroshio current.
We employ this tool on velocity data sets on the Gulf of Mexico obtained
from HYCOM project. The tool identifies the underlying Lagrangian structures in the oceanic currents. In particular invariant manifolds, hyperbolicand non-hyperbolic flow regions are detected. We study the influence of thesestructures on particle motions on the oil spill.
We acknowledge to CESGA for support with the supercomputer FINISTERRAE. Thanks for support by grants: CSIC-OCEANTECH, I-Math C30104, MICINN-MTM2008-03754 and MTM2008-03840-E, and C. MadridSIMUMAT.

80

References[1] S. Bowen, M.S. Lozier, S. F. Gary, and C.W. Bning. Nature 459, 243247, 2009.
[2] K. Katija, and J.O. Dabiri. Nature 460, 624-626, 2009.
[3] Q. Schiermeier. Nature 447, 522-524, 2007.
[4] J. A. Jiménez Madrid, A.M. Mancho. Distinguished trajectories in timedependent vector fields. Chaos 19 (2009), 013111-1-013111-18.
[5] C. Mendoza, A.M. Mancho. The hidden geometry of ocean flows. Physical Review Letters 105 (2010), 3, 038501-1-038501-4.
[6] A.M. Mancho, S. Wiggins, A. Turiel, E. Hernández-Garca, C. López, E.
Garcı́a-Ladona. Nonlinear Processes in Oceanic and Atmospheric Flows.
Nonlinear Proc. Geoph 17 (2010), 3, 283-285.
[7] C. Mendoza, A.M. Mancho, Marie-Hélene Rio. The turnstile mechanism
across the Kuroshio current: analysis of dynamics in altimeter velocity
fields. Nonlinear Proc. Geoph 17 (2010), 2, 103-111.
[8] C. Mendoza, A.M. Mancho. The Lagrangian description of aperiodic
flows: a case study of the Kuroshio Current. preprint (2011).

81

Mathematical studies on the self-motion of
surfactant scrapings at the air-water interface
Keita Iida1),∗ , Hiroyuki Kitahata2,5) ,Satoshi Nakata3) and Masaharu Nagayama4,6)1) Graduate School of Natural Science and Technology, Kanazawa
University,2) Graduate School of Science, Chiba University,3) Graduate School of Science, Hiroshima University,4) Institute of Science and Engineering, Kanazawa University,5) PRESTO, Japan Science and Technology Agency,
6) CREST, Japan Science and Technology Agency∗ e-mail: k.iida@stu.kanazawa-u.ac.jp
Abstract
Several types of autonomous motors are known to exhibit various manners of “self-motion” at air/water interfaces. Spontaneous motion of surfactant scrapings at the water surface (henceforth surfactant-water system)
is a representative example of self-motion, which involves a variety of rich
topics and which was introduced in 1997 [1, 2]. As presented in [2], camphor
scrapings (diameter, ca. 1 mm) exhibit diﬀering manners of sustained motion
on water under almost isothermal and non-equilibrium conditions, of which
the working mechanisms are similar to living organisms. Hence, from thephysicochemical point of view, they have been considered as a novel eﬃcientchemo-mechanical transducer which mimics living things.
From the analytical point of view, we would like to discover the mathematics and physics that underlie the surfactant-water system and to understand the essential mechanisms of the self-motion.* 1) One proposition ona camphor-water system which has been veriﬁed to be true, by means ofexperiments and analyses, is that a pitch-fork bifurcation of travelling pulse
solution occurs depending on the viscosity of the aqueous phase in a onedimensional problem [3]. For the one-dimensional problem, we can now see
such mathematical structure at the point of demarcation between self-motion
and the rest state.
How about the two-dimensional problem? The next step we discuss below is to consider both translational and rotational motion induced by the* 1)

Our real goal is to understand the principles and structures that living and nonliving
systems have in common and to realize the autonomic motor system in living things attheir most fundamental level.

82anisotropic shape of the camphor scraping. Although a mathematical model
consisting of translational and rotational equations, coupled with a reactiondiﬀusion equation, has been proposed in 2005 [4], the characteristics of thetwo-dimensional camphor motion depending on its shape has not yet beenclariﬁed. For example, let us consider a camphor scraping with ellipticalshape. Which direction does the camphor move, in the long axis or in the
short axis? Or, does the elliptical camphor rotate? In fact, we found itdiﬃcult to investigate the above question by brief experiments and simulations due to diﬃculties originating from the elliptical shape (e.g., we have tocompute the elliptic integral in the model).
To tackle the above questions, we apply bifurcation theory and conduct
proper experiments to determine the mode of motion of elliptical camphor
scrapings. In the presentation, we will introduce our mathematical ﬁndings,
numerical simulations and experiments of the camphor-water system.

References[1] A. Mikhailov and D. Meinköhn, Self-Motion in Physico-Chemical Systems Far from Thermal Equilibrium, Springer, Berlin Heidelberg, 484(1997).
[2] S. Nakata and Y. Iguchi, S. Ose, M. Kuboyama, T. Ishii and K.
Yoshikawa, Self-Rotation of a Camphor Scraping on Water: New Insight into the Old Problem, Langmuir, 13 (1997), 4454–4458.
[3] M. Nagayama, S. Nakata, Y. Doi and Y. Hayashima, A theoretical andexperimental study on the unidirectional motion of a camphor disk,
Physica D (Amsterdam), 194 (2004), 151–165.
[4] H. Kitahata and K. Yoshikawa, Chemo-mechanical energy transduction through interfacial instability, Physica D (Amsterdam), 205 (2005),
283–291.

83

A curious bifurcation diagram arising in a scalar
reaction-diffusion equation with an integral constraint
Hisashi OKAMOTO
Research Institute for Mathematical Sciences,
Kyoto University, Kyoto 606-8502* e-mail: okamoto@kurims.kyoto-u.ac.jp
Abstract
We consider uxx + Ru2 = constant + f (x), together with ∫-ππu(x) dx = 0. A

strange bifurcation diagram was found. The external force f is chosen so that u =cos kx becomes a solution (k = 2, 3, 4). Solving the equation by the spectral method,
we obtained non-trivial solutions. In the figure, where k = 2, the solution branchcomes back to the trivial solution at the points indicated by the arrows. As R
increases, the profile of the solution becomes complicated, but they also lose theircomplexity when they come back to the trivial solution.
In my talk, I will explain these phenomena, emphasizing the contrast with thebifurcating solutions appearing in the Kolmogorov flows of 2D Navier-Stokes
equations [1].

References[1] Sun-Chul Kim & Hisashi Okamoto, Vortices of large scale appearing in the
2D stationary Navier-Stokes equations at large Reynolds numbers, Japan J.
Indust. Appl. Math., vol. 27 (2010), 47--71.
84

Application of topological computation
method for global dynamics of an associative
memory model of the Hopfield type
Hiroshi Kokubu1),∗ ,1) Department of Mathematics, Kyoto University / JST CREST∗ e-mail: kokubu@math.kyoto-u.ac.jp
Abstract
A new rigorous topological computation method for global analysis ofdynamical systems was introduced in [1]. I would like to revisit the classical associative memory model of the Hopfield type, and study its global
phase space structure from the point of view of this topological computation
approach.

References[1] Z. Arai, W. Kalies, H. Kokubu, K. Mischaikow, H. Oka, and P. Pilarczyk: A database schema for the analysis of global dynamics of multiparameter systems, SIAM Journal on Applied Dynamical Systems, 8(2009), 757–789.

85

Rigorous verification of equilibria forevolutionary equations- Existence, Uniqueness and Hyperbolicity Kaname Matsue1),∗ ,1) Department of Mathematics, Kyoto University∗ e-mail: k-matsue@math.kyoto-u.ac.jp
Abstract
The aim of this talk is to provide a method for investigating a precise
structure of equilibria of (infinite dimensional) dynamical systems, such as
the local uniqueness and the hyperbolicity.
Recently several rigorous computation methods for verifying stationary
solutions of partial differential equations are studied (e.g. [8] or [10]). Existence and local uniqueness are mainly investigated therein. Moreover, in[2], one of applications of [10], the global structure of all bounded global
solutions for a partial differential equations is studied by using computer assisted proof and an algebraic-topological invariant of invariant sets, called
the Conely index. This invariant gives us the information of an invariant set,
like the existence and its stability.
Hyperbolicity, which is related to the structural stability is known as oneof the most important properties of dynamical systems, and is stronger than
the robustness provided by the Conley index. An equilibrium u of a PDE is
hyperbolic if the linearized operator at u has no spectrum on the imaginary
axis. This property is important from the viewpoint of dynamical systems,
such as the verification of local or global bifurcations (e.g. [3]).
In order to verify the hyperbolicity, one may need analytic methods forinvestigating the spectrum of a linear operator in general. It is sometimes
the case, however, that a linearized operator is self-adjoint, in which case theinvertibility of the linearized operator may be sufficient to yield the hyperbolicity of an equilibrium. Even if the linearized operator is not self-adjoint, we
can prove the hyperbolicity of equilibria for a broad class of infinite dimensional dynamical systems, called dynamical systems generated by parabolic
evolutionary equations ([5]).
The author shows that the specific Conely index and the suitable hypothesis of the linearized operator yield the local uniqueness and the hyperbolicity
of an equilibrium. Using the preceding work by Nakao and his collaborators([7]) and by Zgliczyński and Mischaikow ([10]) or the author ([6]), we can86verify the hyperbolicity of equilibria of parabolic evolutionary equations bycomputer assisted analysis.
At the end of this talk, the author will show several sample rigorous
numerical results about the existence, local uniqueness and hyperbolicity ofequilibria.

References[1] C.Conley, Isolated Invariant Sets and the Morse Index, CBMS Reg. Con.
Ser. Math., Vol.38, Amer. Math. Soc. Providence, RI, 1978.
[2] S. Day, Y. Hiraoka, K. Mischaikow and T. Ogawa, Rigorous Numerics for
Dynamics: A Study of the Swift-Hohenberg Equation, SIAM J. Applied
Dynamical Systems, 4(2005), 1–31.
[3] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, 1983.
[4] J.K. Hale, Asymptotic behavior of dissipative systems, Mathematical
Surveys and Monographs, 25. AMS, Providence, 1988.
[5] D. Henry, Geometric Theory of Semilinear Parabolic Equations,
Springer-Verlag, 1981.
[6] K. Matsue, Rigorous verification of equilibria for evolutionary equations - Existence, Uniqueness and Hyperbolicity - Part 1. Gradient case,preprint.
[7] M.T. Nakao, K. Hashimoto and Y. Watanabe, A Numerical Method to
Verify the Invertibility of Linear Elliptic Operators with Applications to
Nonlinear Problems, Computing, 75(2005), 1–14.
[8] M.T. Nakao and N. Yamamoto, Numerical verifications for solutions
to elliptic equations using residual iterations with higher order finite
element, J. Comput. Appl. Math., 60(1995), 271–279.
[9] K.P. Rybakowski, The Homotopy Index and Partial Differential Equations, Springer-Verlag, 1987.
[10] P. Zgliczyński and K. Mischaikow, Rigorous numerics for partial differential equations: The Kuramoto-Sivashinsky Equation, Found. Comput.
Math. 1(2001), 255–288.

87

Role of unstable symmetric solutions insymmetry restoring process
Tsuyoshi Mizuguchi1),2)∗ , Makoto Yomosa1)
1) Department of Mathematical Sciences, Osaka Prefecture University.
2) PRESTO, Japan Science and Technology Agency (JST).
∗ e-mail: gutchi@ms.osakafu-u.ac.jp
Abstract
In the study of dynamical systems, symmetries are one of the most important concepts for characterising and/or classifying systems and solutions.
For example, some bifurcations involve spontaneous breakdown or restoration of symmetry. Several chaotic systems having a discrete symmetry often
exhibit a bifurcation that breaks or restores the symmetry while maintaining
its aperiodic behaviour, i.e., an attractor merging crisis (AMC).
On the asymmetric side of the AMC bifurcation point, a pair of strange
attractors exists, one of which is chosen by an initial condition. At thebifurcation point, they merge with each other. On the symmetric side, asingle strange attractor that satisﬁes the symmetry is observed. Let us focus
on singularities near the AMC bifurcation point. Intermittent behaviours
are reported in several systems. Namely, the inter-burst interval time tends
to diverge when the system approaches the bifurcation point. If a state istranslated to a motion in space, the singularity corresponds to divergenceof the diﬀusion coeﬃcient of the motion. Therefore, one can forecast thebifurcation by observing these quantities.
The singularity described above is observed when the system approaches
from the symmetric side. On the other hand, from the asymmetric side,a speciﬁc unstable symmetric solution so-called mediating unstable periodic
orbit, or simply a mediating solution M plays an important role to characterize the symmetry restoration process. The mediating solution is deﬁned
as that with which a pair of strange attractors simultaneously collides at the
AMC point when the system approaches from the asymmetric side. Near thisAMC point, the outermost part of the strange attractor in each parameter
passes close to solution M . Therefore, the distance dM between the strange
attractor and M is an appropriate quantity to characterise the closeness tothe bifurcation point.
This approach is also characterised by a duration in which a state isclose to solution M . We deﬁne a state whose distance to M is less than athreshold value as a temporarily sticking state, and we focus on the maximum
88

duration τM of this temporarily sticking state during the measurement time.
A logarithmic singular behaviour of τM with periodic modulation is predicted
near the bifurcation point.
In this paper, we exemplify several dynamical systems exhibiting AMC
transition, such as Ginzburg Landau type equation with external force termor coupled Stuart Landau type equation. The distance dM and the maximum duration τM are directly measured near the bifurcation point. A linear
decrease in dM and the logarithmic divergence of τM are conﬁrmed numerically. We suggest a time series analysis method to measure these quantities
without knowing the precise bifurcation point or concrete function form of
the mediating solution M .

References[1] M. U. Kobayashi, T. Mizuguchi, “Chaotically oscillating interfaces in aparametrically forced system”, Phys. Rev. E, 73 (2006) 016212.
[2] M. U. Kobayashi, T. Mizuguchi, “Chaotic interfaces in a parametricallyforces system”, Prog. Theor. Phys. Supplement, 161 (2006) 228–231.
[3] Y. Morita, N. Fujiwara, M. U. Kobayashi and T. Mizuguchi, “Scytale
decodes chaos: A method for estimating unstable symmetric solutions”,
Chaos, 20 (2010) 013126.

89

A Diffusion Equation Representing the Dynamics of the
Jakarta Composite Index (JSX)1),2)*Edi Cahyono1), 3), Buyung Sarita

1) Graduate Program, Universitas Haluoleo2) Department of Mathematics FMIPA, Universitas Haluoleo
3) Faculty of Economics, Universitas Haluoleo* e-mail: edi_cahyono@innov-center.org
Abstract
We cosider the dynamics of the Jakarta Composite Index (JSX). The dynamics
is presented as time series. The dynamics may be charaterized by the trend and theprobability density function (pdf) based on the assumption that at a fixed time it isnormally distributed with respect to ‘spatial’ variable. Details of the trend and the
pdf for the case of the dynamics of JSX will be presented in the succesive article[2]. Following the work of [1], we seek a type of diffusion equation which has
fundamental solution representing the trend and the probability density function ofsuch dynamics.

References[1] E. Cahyono, J.R. Juliana, R. Raya: Asset value dynamics and thefundamental solution of a modified heat equation, 4th Intl Conf. Math. Stat.,
Bandar Lampung, Indonesia (2009), 312--318.
[2] B. Sarita, F. Arisanti, E. Cahyono: The trend and probability density functionof Indonesia Stock Index (IHSG), submitted.

90

The Trend and Probability Density Function of the
Jakarta Composite Index (JSX)1),2)*

Buyung Sarita

3)1), 3), Femy Puspita Arisanty , Edi Cahyono

1) Graduate Program, Universitas Haluoleo
2) Faculty of Economics, Universitas Haluoleo3) Department of Mathematics FMIPA, Universitas Haluoleo* e-mail: buyungsarita_buyung@yahoo.com
Abstract
In this paper we discuss the dynamics of the Jakart Composite Index (JSX),
which is related to the work of [3]. The dynamics is presented as time series. To
predict the dynamics is still difficult. In general it is almost impossible to predict
such dynamics for the case of high frequency data [2]. Hence, we do not predict the
dynamics. Rather, we seek the trend and the probability density function based on
the assumption that the dynamics is normally distributed with respect to time
variable. According to [1], the trend and the probability density function of such
dynamics is the fundamental solution of a type of diffusion equation.

References[1] E. Cahyono, J.R. Juliana, R. Raya: Asset value dynamics and thefundamental solution of a modified heat equation, Proc. 4th Intl Conf. Math.
& Stat., Bandar Lampung, Indonesia (2009), 312--318.
[2] T. Ito, K. Sato: Exchange rate change and inflation in post-crisis Asian
economies: Vector autoregression analysis of exchange rate pass-through, J.
Money, Credit and Banking, (2008) 40(7), 1407--1437.
[3] B. Sarita, Agusrawati, L.D. Suryadi, S. Wahab, E. Cahyono: Dynamics ofgold price relative to Indonesian rupiah, Proc. 3rd Intl. Conf. Quantitative
Methods, Bandar Lampung, Indonesia (2010) 42--46.

91                          !"!! !#$    " %      &&    ' "  "     % && ( &&&          " " & ! )"  " " & &&     % "    *   "    "%    !
$ "           "   "  " "         "    " &   ! "  " "  &           "        " " "        % "         + ,!    " "           -   "      "     +.,! /-   " %& *   "            
 " 0   "&    +1,!   
" %& *                 "  " "  &          
 "      "    +2,! )"   "   "     "%( & &&   !
½µ    + , !   )! /%  )! 3& 4        "         "  -    " 5   "  6& ! 78 !2792:9 (.;;9+., !    ! <" -    "      "        <" ! 9
;21 . (.;;9+1, )! /% !    /! =  $   &&&  %  "  -    " >  &! 71 !12 17; (.;92

+2, !    )! /%         "         ?$4 5   4 "  ! :  ! @9 :9 (.;93

Truncation error analysis for particle methods
Masato Kimura1),∗ , Kiyohiro Ishijima2)
1) Institute of Mathematics for Industry, Kyushu University
2) Graduate School of Mathematics, Kyushu University∗ e-mail: masato@imi.kyushu-u.ac.jp
Abstract
Meshfree particle methods such as moving particle semi-implicit (MPS)
method and smoothed particle hydrodynamics (SPH) method have beenmuch developed in recent years and established as powerful tools in the fields
of science and engineering, since they require no troublesome mesh generation and have flexible adaptation to changing shape of domains. On theother hand, there are few mathematical error analysis on their accuracy as
discretization methods for partial differential equations.
In [1], we proposed a new quantity concerning the distribution of theparticles and gave a truncation error estimate for a discrete gradient operatorof meshfree particle methods with weight functions, particularly one of MPS
method.
Let Ω be a subset of Rm (m ≥ 2) with a piecewise smooth boundary. We
suppose thatxi 6= xj (i 6= j)
X = {xi }N
i=1 ⊂ Ω,
is given, and we call xi a particle. In order to measure the uniformness of the
distribution of X in Ω, we introduce the following quantities. If a divisionω = {ωi }N
i=1 of Ω satisfies the condition:
ωi ⊂ Ω,

ωi : open,

ωi ∩ ωj = ∅ (i 6= j),|ωi | =

1|Ω|N

we call ω N -equivolume partition of Ω.
For given N particles X and an N -equivolume partition ω in Ω, we define
the equivolume partition radius δ∞ (X) by
δ∞ (X) := inf δ∞ (X, ω),
ω

(1)

where
δ∞ (X, ω) := max sup |x − xj |,1≤j≤N x∈ωjand the infimum in (1) is took over all N -equivolume partition of Ω. The
equivolume partition radius δ∞ (X) is a quantity which stands for a kind of94

distance between particles. If X distributes in Ω almost uniformly, there
exist C1 > 0 and C2 > 0 independent of N such that
1

1

C1 N − m ≤ δ∞ (X) ≤ C2 N − mholds.
In MPS method, we use a weight function
1− 1, 0(r > 0),wMPS (r) = C maxr
R
where C > 0 is chosen as Rm w(|x|) dx = 1. Moreover, for a > 0, we
define wa (r) := a−m w( ar ). The support of wa is the ball Ba (xi ) := {x ∈Rm ; |x − xi | ≤ a}. We call a the influence radius. Using the scaled weightfunction wa , we define the particle density di (X, a) and the particle density
deviation qi (X, a) as follows:|Ω| X adi (X, a) :=w (|xj − xi |),qi (X, a) := di (X, a) − 1.
(2)
N j6=iFor a function φ defined on X with φi = φ(xi ) (i = 1, 2, · · · , N ), thefollowing discrete gradient operator is used in MPS method [2]:m|Ω| X (φj − φi )ah∇φii :=
2 (xj − xi ) w (|xj − xi |)N j6=i |xj − xi |Even for a function u(x) defined on Ω, we define h∇ui similarly with ui =u(xi ). Our main result for m = 3 is as follows.
Theorem 1. Let m = 3 and let X = {xj }N
j=1 be a set of particles in Ω. If
xi ∈ X and a > 0 satisfy the condition Ba (xi ) ⊂ Ω, for u ∈ C 3 (Ω), thefollowing inequality holds:δ∞ (X)
2|h∇uii − ∇u(xi )| ≤ C1
+ a + C2 qi (X, a).
aReferences[1] K. Ishijima, M. Kimura: Truncation error analysis of finite difference
formulae in meshfree particle methods, Trans. Japan Soc. Indust. Appl.
Math, 20 (2010), 165-182. (in Japanese)[2] S. Koshizuka, Y. Oka: Moving particle semi-implicit method for fragmentation of incompressible fluid, Nuclear Science and Engineering, 123(1996), 421434.
95

A numerical method for multiphasevolume-preserving mean curvature flow
Elliott Ginder1),∗ , Seiro Omata2) , Karel Švadlenka3)
1) Grad. Sch. of Nat. Sci. and Tech., Kanazawa University
2) Institute of Science and Engineering, Kanazawa University
3) Institute of Science and Engineering, Kanazawa University∗ e-mail: eginder@polaris.s.kanazawa-u.ac.jp
Abstract
We introduce a method for computing multiphase motion by volumepreserving mean curvature ﬂow. The method is variational and can be implemented in any dimension, with interfaces in the two dimensional setting
described by closed curves, each of which encloses a prescribed area. The
corresponding physical phenomena for curvature driven motions are numerous and include the grain boundary growth of crystals and ternary alloys,
and applications also arise in image processing and biology.
Frequently occurring in practical applications, curvature driven interfacial
motions are a recurring subject of study. As such, analytical and numericalmethods have been proposed for their treatment. Standard approaches include front tracking [8], level set [5], and phase ﬁeld [1, 6] methods, but thesetechniques are limited to two phases whenever they are extended to handle
the case of volume preservation (see, [2, 3]).
Our approach to realizing multiphase volume-preserving mean curvature
ﬂow is based on an idea of Merriman, Bence and Osher (MBO)[4]. The MBO
algorithm approximates mean curvature ﬂow by diﬀusing, for a small time,
the characteristic function of the set enclosed by the interface. In particular,
the characteristic function of the set is taken as the initial condition forthe heat equation, which is then solved for a short time under Neumannboundary conditions. The location of the evolved interface is then given by
the 1/2 level set of the solution. The initial condition for the heat equationis renewed by truncation (that is, by taking the characteristic function of the
1/2 level set), and the process is repeated.
The main beneﬁt of this approach is that the diﬀusion process is able
to automatically treat the occurrence of singularities, and there is no needdirectly compute curvatures. This method is given for the multiphase case in[4], and for a 2-phase volume-preserving method, see [7]. However, the ideain [7] does not apply to the multiphase case, and so we reformulate the ideaof [4], and extend it to the case of multiphase motion by volume-preserving
mean curvature ﬂow.
96

The prescribed phase volumes are maintained by use of a variational technique known as the discrete Morse flow. Time is discretized in this method
and a functional is introduced for each time step. The functionals correspondto the discrete Morse ﬂow for a vector valued heat equation, with additional
penalty terms for the prescribed phase volumes. Minimizing the functionals
decreases their penalties, so that each phase maintains its prescribed volume.
Similar to the MBO algorithm, we then employ a truncation to obtain theevolution of the interfaces, and we repeat.
We will explain the details of our algorithm and show results of its application in the two dimensional setting.

References[1] L. Bronsard, F. Reitich: On three-phase boundary motion and the singular limit of a vector-valued Ginzburg-Landau equation, Arch. Rational
Mech. Anal., 124 (1993), pp. 355-379.
[2] L. Bronsard, B. Soth: Volume preserving mean curvature ﬂow as a limitof a nonlocal Ginzburg-Landau equation, SIAM J. Math. Anal., 28 No.
4 (1997), pp. 769-807.
[3] S. Esedoglu, P. Smereka: A variational formulation for a level set representation of multiphase ﬂow and area preserving curvature ﬂow, Communications in Mathematical Sciences, 6 No. 1 (2008), pp. 125-148.
[4] B. Merriman, J. Bence, S. Osher: Motion of multiple junctions: a levelset approach, J. Comp. Phys., 112 No. 2 (1994), pp. 334-363.
[5] S. Osher, J. A. Sethian: Fronts propagating with curvature-dependent
speed: Algorithms based on Hamilton-Jacobi formulations, J. Comp.
Phys., 79 (1988), pp. 12-49.
[6] J. Rubinstein, P. Sternberg: Nonlocal reaction-diﬀusion equations andnucleation, IMA J. Appl. Math., 48 (1992), pp. 248-264.
[7] S. J. Ruuth, B T. R. Wetton: A simple scheme for volume-preserving
motion by mean curvature, J. Scientiﬁc Computing, 19 (2003), pp. 373384.
[8] J. E. Taylor: The motion of multiple-phase junctions under prescribedphase-boundary velocities, J. Diﬀerential Equations, 119 (1995), pp.
109-136.

97

A High-Order Discontinuous Galerkin
Method for Elliptic Interface Problems
Min-Hung Chen1),∗ , Rong-Jhao Wu2)
1) Department of Mathematics, National Cheng Kung University, Taiwan2) Department of Aeronautics and Astronautics, National Cheng Kung
University, Taiwan∗ e-mail: minhung@gmail.com
Abstract
In this work, we implement the high-order local discontinuous Galerkinmethod [1] with the numerical ﬂux proposed by Guyomarc’h et. al. [3], tosolve the elliptic interface problem.
Instead of the triangular elements and the P k -polynomial space used in[3], we use the quadrilateral elements and the Qk -polynomial space. To treat
the curved interface and boundaries, we use the transﬁnite blending mappings[2] to construct curved elements. The discretization generates a symmetricsystem and the system can be solved directly with standard methods or
reduced into smaller systems with matrix reduction techniques.
The numerical experiments show the h and p convergence properties ofthe our high order scheme. Moreover, our numerical experiments show thatthe high-order method is more eﬃcient than the low-order one. To achieve
the same magnitude of accuracy, the degree of freedom of the generated
system using the coarsest mesh is much smaller than the one using the ﬁnest
mesh. As a result the CPU time required to solve the system using thecoarsest mesh is much smaller, too.

References[1] Cockburn, B.; Shu, C.-W.: The local discontinuous Galerkin method
for time-dependent convection-diﬀusion systems. SIAM J. Numer. Anal.
vol. 35, no. 6, pp. 2440-2463 (1998).
[2] Gordon, W. J.; Hall, C. A.: Transﬁnite element methods: blendingfunction interpolation over arbitrary curved element domains. Numer.
Math. vol. 21, no. 2, pp. 109-129 (1972).
[3] Guyomarc’h, G.; Lee, C.O.; Jeon, K.: A discontinuous Galerkin method
for elliptic interface problems with application to electroporation, Commun. Numer. Methods Engrg. (2008).

98

A Coarse-Grain Parallel Scheme for Solving Poisson
Equation by Chebyshev Pseudospectral Method
Teng-Yao Kuo1); , Hsin-Chu Chen2) , Tzyy-Leng Horng3)
1) Ph. D. Program in Mechanical and Aeronautical Engineering, Feng Chia
University, Taichung, Taiwan2) Department of Computer and Information Science, Clark Atlanta
University, Atlanta, GA 30314, USA
3) Department of Applied Mathematics, Feng Chia University, Taichung,Taiwane-mail: troyguo@yahoo.com.tw
Abstract
Poisson equation is frequently encountered in mathematical modeling forscienti…c and engineering applications. Fast Poisson numerical solvers for 2D
and 3D problems are, thus, highly requested for their simulations. In this paper, we consider solving the Poisson equation r2 u = f (x; y) in the Cartesian
domain = [ 1; 1] [ 1; 1], subject to mixed Dirichlet / Neumann boundary conditions, discretized with the Chebyshev pseudospectral method. The
main purpose of this paper is to propose a re‡exive decomposition scheme
for orthogonally decoupling the linear system obtained from the discretization into independent subsystems via the exploration of a special re‡exivity
property inherent in the second-order Chebyshev di¤erentiation matrix. The
decomposition not only yields more e¢ cient algorithm but introduces coarsegrain parallelism. This approach can be applied to more general problems
and to other partial di¤erential equations discretized with the Chebyshev
pseudospectral method as well, so long as the discretized problems possessre‡exive symmetry. Numerical examples and the performance of numericalexperiments are presented to demonstrate the validity and advantage of the
proposed approach.
Keywords Poisson equation, Chebyshev pseudospectral method, Chebyshev di¤erentiation matrix, Coarse-grain parallelism, Re‡exivity property

References[1] A. L. Andrew: Solution of equations involving centrosymmetric matrices, Technometrics, 15(2) (1973), 405–407.
[2] A. L. Andrew: Eigenvectors of certain matrices, Linear Algebra Appl.,
7 (1973), 151–162.
99

[3] A. Cantoni, P. Butler: Eigenvalues and eigenvectors of symmetric centrosymmetric matrices, Linear Algebra Appl., 13 (1976), 275–288.
[4] C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zang: Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1987.
[5] H-C. Chen: Increasing parallelism in the …nite strip formulation: static
analysis, Neural, Parallel & Scienti…c Computations, 2 (1994), 273–298.
[6] H-C. Chen, A. Sameh: A matrix decomposition method for orthotropic
elasticity problems, SIAM J. Matrix Anal. Appl., 10(1) (1989), 39–64.
[7] H-C. Chen, T-L Horng, Y-H. Yang: Re‡exive decompositions for solving
Poisson equation by Chebyshev pseudospectral method, Proceedings of
Neural, Parallel, and Scienti…c computations, 4 (2010), 98–103.
[8] R. Peyret: Spectral methods for incompressible viscous ‡ow, Applied
Mathematical Sciences 148, Springer-Verlag, 2002.
[9] H. Dang-Vu, C. Delcarte: An accurate solution of the Poisson equationby the Chebyshev collocation method, J. Comput. Phys., 104 (1993),
211–220.
[10] U. Ehrenstein, R. Peyret: A Chebyshev collocation method for NavierStokes equations with applications to double-di¤usive convection, Int’lJ. for Numerical Methods in Fluids, 9 (1989), 427–452.
[11] G.H. Golub, C.F. Van Loan: Matrix Computations, The Johns Hopkins
University Press, Maryland, USA, 1983.
[12] D. B. Haidvogel: The solution of Poisson equation by expansion in
Chebyshev polynomials, J. Comput. Phys., 30 (1979), 167–180.
[13] J.d.J. Martinez, P.d.T.T. Esperanca: A Chebyshev collocation spectral
method for numerical simulation of incompressible ‡ow problems, J. of
the Braz. Soc. of mech. Sci. & Eng., XXIX(3) (2007), 317–328.
[14] Lloyd N. Trefethen: Spectral Methods in MATLAB, SIAM, Philadelphia, 2000.
[15] J.R. Weaver: Centrosymmetric (cross-symmetric) matrices, their basic
properties, eigenvalues, and eigenvectors, The American Mathematical
Monthly, 92(10) (1985), 711–717.

100

On the stability of an oscillation-free ADI
method in laser beam propagation
computations
Qin Sheng1) , Hai-Wei Sun2),∗1) Department of Mathematics, Center for Astrophysics, Space Physics andEngineering Research, Baylor University, USA
2) Department of Mathematics, University of Macau, Macao, China∗ e-mail: hsun@umac.mo
Abstract
If a slowly varying envelope approximation of the light beam, such as alaser, is concerned, then one of the most effective mathematical models used
can be the paraxial Helmholtz equation,2iκuz = uxx + uyy + g(u),∀(x, y, z) ∈ D,

(1)√
where i ≡ −1 is the imaginary unit, κ is the wave number, z is the beam
propagation direction, x, y are dimensional directions perpendicular to thelight, u is the complex envelope of the wave to be investigated, g is the source,
and D = {a < x < b, c < y < d, z > z0 } . Equation (1) is often called theone-way Helmholtz equation and can be viewed as a dimensionless simplification of Maxwell’s field equations. Moreover, (1) has been widely usedin laser optical material researches, particularly in crystal or liquid crystal
photorefractive and negative indexed (or left-handed) material developments.
Since κ can be extremely large in optical applications, the complex envelope is highly oscillatory. Consequently, higher efficiencies in solution computations for (1) are difficult to achieve, since mesh steps cannot be unrealistically small. Although for the spectrum, Gaussian beam methods andintegral transformations may possess certain advantages in the situation, themain challenge pertains in balancing the algorithmic simplicity and accuracy.
This motivates a recent study of the eikonal transformation based numerical methods. For simplicity in the following discussion, we remove thesource g(u) from (1). Assume φ(x, y, z) and ψ(x, y, z) to be the amplitude
and phase functions of the underlying highly oscillatory wave respectively.
Based on the theory of geometrical optics, we claim thatu(x, y, z) = φ(x, y, z)ei·κψ(x,y,z) ,101∀(x, y, z) ∈ D.

(2)

By substituting (2) into (1), we immediately acquire the following coupledeikonal differential equations in D, i.e.,φz = α (ψxx + ψyy ) + f1 ,
ψz = β (φxx + φyy ) + f2 ,

(3)
(4)

where
α=

i
φ
1
1h
(ψx )2 + (ψy )2 .
, β = − 2 , f1 = φx ψx + φy ψy , f2 =
2
2κ φ
2

(5)

Note that φ 6= 0, solutions of (3) and (4) are real and non-oscillatory
even when κ is extremely large. Therefore, mesh steps anticipated for solving
them need not to be particularly small. Needless to say, this implies a terrificpotential for fast and highly efficient computations of the numerical solutionof (1). Discussions of continuing discretizations of (3) and (4) can be foundin numerous recent publications, in particular, with important electro-opticaland material engineering applications.
Although various kinds of finite difference schemes based on (3) and (4)have been accomplished and used in practical laser computations, their numerical stability has been a serious concern in many applications. In fact,
the algorithms comprised are often unstable, especially when a large beam
propagation distance is considered. This proposes an extremely importantquestion on the applicability of the eikonal transformation based schemes inmodern laser beam computations. Therefore, it becomes more imminently
for us to study the stability of the popular oscillation-free split step method.
In this talk, we study the stability of an ADI finite difference method
for solving highly oscillatory paraxial Helmholtz equations in laser optics
applications. We show that the ADI method is asymptotically stable whenwave numbers anticipated are large.

102

Flows in thoracic aorta with torsion
Hiroshi Suito1),4)∗ , Takuya Ueda2),4) , Daniel Sze3),
1) Graduate School of Environmental Science, Okayama University, Japan
2) Department of Radiology, Saint Luke’s Hospital, Japan
3) Department of Radiology, Stanford University School of Medicine, USA
4) CREST, JST∗ e-mail: suito@ems.okayama-u.ac.jp
Abstract
Thoracic endovascular aortic repair (TEVAR), or stent-graft treatment,
has become widely accepted as an important option for treatment of thoracic
aortic diseases. Many studies have proven the safety and eﬃcacy of TEVAR
with satisfactory short-term to mid-term outcomes. Nevertheless, even if theinitial TEVAR treatment technically succeeds, some patients show recurrence and progression of disease many years after treatment[5, 6]. Based onlong-term follow-up examinations, such long-term morphological change andeﬀects of hemodynamic ﬂow apparently interact synergically. Constant pulsatile hemodynamic eﬀects from blood ﬂow apparently induce degenerationof the underlying aorta to cause its morphological change and to induce minor morphological changes that alter the hemodynamic state. These changes
ultimately engender long-term adverse events. For this study, which investigates the eﬀect of vascular hemodynamics on long-term adverse events, weconsidered patient-speciﬁc models of the thoracic aorta as constructed fromCT scans.
As a procedure for the reconstruction of aorta morphology from medical
images, median axis transform technique[4] is applied to extract the (x, y, z)coordinate of the centerline and radius at position s from CT images, wheres is a length along the centerline from the proximal end of aorta. Then aﬁnite diﬀerence mesh system (ξ, η, ζ) is generated, one coordinate axis ofwhich is nearly parallel to the aorta centerline. The ζ−axis is set not to bestrictly parallel but to be nearly parallel to the centerline. This treatment
avoids numerical instabilities arising from severe skewness of ﬁnite diﬀerence
meshes. On each mesh point, the characteristic function is computed, whichtakes the value 0 in the aorta and value 1 out of the aorta. It varies graduallynear the boundary with certain smoothness.
As governing equations, incompressible viscous Navier-Stokes equations
are used with the continuity equation. The governing equations are discretized on the ﬁnite diﬀerence mesh and SMAC method is applied for time103

advances. At the proximal end, a pulsating velocity proﬁle is given, whereas
the constant pressure is given at the distal end.
It is apparent from numerical results that the stress distributions arestrongly dependent to the aorta morphology. From comparison of the ﬂow
ﬁelds of several examples of thoracic aorta morphologies, it has been inferredthat torsions of the centerline are a candidate of the critical parameter of theexistence of swirling ﬂow in diastole phases. In order to examine the eﬀect of
torsion on the ﬂows in curved tubes, numerical tests in simpliﬁed geometries
are performed.
Through computations of the ﬂows in realistic aorta morphologies andin simpliﬁed spiral tubes, the importance of the torsion of aorta has beenshown. The values of torsion can vary widely among patients. Our results
suggest that the torsion is a more important parameter in considering therelation between aorta morphologies and the hemodynamics in them.

References[1] C.Y. Wang: On the low-Reynolds-number ﬂow in a helical pipe, J. Fluid
Mech., 108 (1981), 185–194.
[2] M. Germano: On the eﬀect of torsion on a helical pipe ﬂow, J. Fluid
Mech., 125 (1982), 1–8.
[3] H. Fujita, H. Kawahara, and H. Kawarada: Distribution theoretic approach to ﬁctitious domain method for Neumann problems, East-West
J. Numer. Math., 3, 2 (1995), 111–126.
[4] G.D. Rubin, D.S. Paik, P.C. Johnson, and S. Napel: Measurement of theaorta and its branches with helical CT, Radiology, 206 (1998), 823–829.
[5] M. Tillich, R.E. Bell, D.S. Paik, D. Fleischmann, M.C. Soﬁlos, L.J.
Logan, and G.D. Rubin: Iliac arterial injuries after endovascular repairof abdominal aortic aneurysms: correlation with iliac curvature anddiameter, Radiology, 219 (2001), 129–136.
[6] T. Ueda, D. Fleischmann, G.D. Rubin, M.D. Dake, and D.Y. Sze: Imaging of the thoracic aorta before and after stent-graft repair of aneurysms
and dissections, Semin. Thorac. Cardiovasc. Surg., 20, 4 (2008), 348–357.
[7] H. Suito, T. Ueda, M. Murakami, and G.D. Rubin: Vortex dynamics inthoracic aortic aneurysms, ECCOMAS CFD 2010, Laboratorio Nacionalde Engenharia Civil, Lisbon, Portugal, (2010).

104

A Pore-scale Network Flow Model for Two Phase Flow: The Genesis
and Migration of Gas Phase
Koukung Alex Chang      
e-mail: chang@mail.npue.edu.tw
Abstract
We present a network flow model to compute the transport, through a pore network, of acompositional fluid consisting of water with a dissolved hydrocarbon gas. Interest in this research is
to develop simulation capability at the pore scale for specific application in order to investigate
reactive flow resulting from CO2 sequestration and capture pore-level flow detail.
The model presented here contains two parts: one describes fluid transport equations which
follow mass conservation law; the other is reactive equations which represent phase equilibrium
conditions. Combining these two parts, this model forms a nonlinear system which contains six
unknowns as well as six equations. The challenge of this model is that the mass conservation
equations and the phase equilibrium conditions be solved simultaneously while the algorithm
maintains enough flexible to allow for the spontaneous creation and extension of the gas phase.
The numerical model employs sequential solution for the pore pressures, saturations, andconcentrations each time-step. The challenge in computing saturations and concentrations is tosatisfy both mass conservation and phase equilibrium. This is done sequentially, updating species
concentrations first due to transport followed by a phase-change/equilibrium (flash) calculation forequilibrium concentration and saturations in each pore.
Computational tests with constant injected flow rate were run on a “quarter 5-spot” flow
scenario to model horizontal flow; the solution obeys the symmetry of the quarter-5-spot pattern.
Under constant temperature conditions, the concentration of CO2 dissolved in the liquid phasethroughout the medium is linearly related to the liquid pressure. In addition, the saturation of the gas
phase throughout the medium is also linearly related to the liquid pressure

1051

Convective transport and stability in films ofbinary mixtures
Santiago Madruga1),∗ , Fathi Bribesh2) Uwe Thiele3)
1) Universidad Politécnica de Madrid, ETSI Aeronáuticos, Plaza Cardenal
Cisneros 3, 28040 Madrid, Spain2), 3) School of Mathematics, Loughborough University, Loughborough,Leicestershire, LE11 3TU, UK∗ e-mail: santiago.madruga@upm.es
Abstract
Thin polymer films are increasingly used in advanced technological applications. The use of these films as coatings is often limited by their lack ofstability due to their wettability properties on the substrates. The instabilities may be employed to create complex morphologies for polymeric functional layers and their development is relatively well understood for single
component liquids. However, in many relevant applications the film consistsof a binary mixture such as a polymer blend. For such systems the dynamics
of the decomposition within the film and of the dewetting of the film itself
may couple. This allows for new pathways of structuring like decomposition
induced dewetting [1] that do not exist in single component films.
We present a model for the evolution of films of isothermal binary liquid
mixtures with a free evolving surface [2]. The model is based on model-H[3] supplemented by appropriate boundary conditions at the free surface andthe solid substrate. The equations account for the coupled transport of theconcentration of a component (convective Cahn-Hilliard equation) and themomentum (Korteweg-Navier-Stokes equation). The inclusion of convective
motion makes surface deflections possible, i.e., the model allows to study
couplings between the decomposition of the mixture and the evolving surface
corrugations.
After determining homogeneous and vertically stratified base states of freesurface films of polymer mixtures we analyse their linear stability with respect
to lateral perturbations [4]. For purely diffusive transport, an increase infilm thickness either exponentially decreases the lateral instability or entirely
stabilizes the film. The inclusion of convective transport leads to a further
destabilization as compared to the purely diffusive case. In some cases theinclusion of convective transport and the related widening of the range ofavailable film configurations (films are then able to change its surface profile)change the stability behavior qualitatively.
106

We show two dominant driving mechanisms for the convective motion
in binary mixtures with diffuse interface: Marangoni driving for energetic
biased surfaces, and Korteweg driving for neutral surfaces.
In addition, we discuss the role of composition for off-critical mixtureson surface deflections, and the dependence of the instability on parameterssuch as the Reynolds number, the surface tension number and the ratio of
velocities of convective and diffusive transport.
S. Madruga acknowledges support by a Marie Curie European Reintegration Grant (PERG04-GA-2008-234384) within the 7th European Community
Framework Programme and by Universidad Politécnica de Madrid throughgrant AL11-PAC-01. U.Thiele acknowledges support by the EU via the ITN
MULTIFLOW (PITN-GA-2008-214919).

References[1] R. Yerushalmi-Rozen, T. Kerle, and J. Klein. Science. 285, 1254–1256(1999).
[2] U. Thiele, S. Madruga, and L. Frastia. Phys. of Fluids. 19, 122106, 2007.
[3] D.M. Anderson, G.B. McFadden, and A.A. Wheeler. Ann. Rev. Fluid
Mech. 30, 139-165 (1998).
[4] S. Madruga and U. Thiele. Phys. of Fluids. 21, 062104, 2009.
[5] S. Madruga and U. Thiele. Eur. Phys. J. Special Topics. 192, 101-108(2011).

107

Longshore Wave Breaker with Reflected Beach
S.R. Pudjaprasetya1),∗ , F.P.H. van Beckum2) , Elis Khatizah3)1) Industrial and Financial Mathematics Research Group,Fac. of Math. and Natural Sciences, Bandung Institute of Technology, Indonesia2) Dept. of Applied Math., University of Twente, the Netherlands
3) Mathematics Dept., Bogor Agricultural University, Indonesia∗ e-mail: sr pudjap@math.itb.ac.id

Abstract
We consider waves pass over a seabed covered with parallel bars as wave
reﬂector. In the previous research [1] an optimal dimension of the submerged
parallel bars is obtained. The more bars installed the more incident wave
amplitude is reduced. In that study we assume the beach is ideal and canabsorb all transmitted wave energy. When there is a reﬂected beach, thesituation is changed because there is a reﬂected wave from the beach. This
reﬂected wave will interact with other waves and will change the long-term
wave interaction behavior. Literatures suggest that submerged bars in frontof beach can either provide shelter or worsen the hazard of incident wave, depending on the phase relation between the bars and the shoreline reﬂection.
In this paper, we study the eﬀectiveness of submerged bars as breakwater
when there is a hard-wall beach on the right. We consider the shallow waterequation in its Riemann invariant form with variables surface elevation anddischarge. Applying the method of characteristic will give us exact formulas
for surface elevation and discharge at any time. In this way we can then obtained the long-term wave interaction behavior. It turns out that the phase
diﬀerence between right running waves will determined whether the wave
will superpose destructively or constructively before they hit the hard-wallbeach. It turns out that distance Lb between submerged bar and the beachplays an important role, because it determines the phase relation between
interacting waves. We obtain formulas for the safest case
√
1 λ1
1
gh1
1 π
= (n + )π,(1)
Lb = (n + ) = (n + )
2 k1
2 2
2
ω
and the most dangerous case
√
π
λ1
gh1
Lb = n = n
= nπ
.
k1
2
ω

(2)

Figure 1 shows results from computation for one-bar submerged bar that
conﬁrm the formulas above. Basically, we can determine the safest distance
Lb for any formation of submerged rectangular bars. Simulations for two-bar
submerged bar and periodic rectangular bar are given.

108

2
1.5
1
0.5
0
−0.5
−1
−1.5
−2

0

20

40

60

80

100

Figure 1: The ampliﬁcation curve of surface elevation at the beach as functionsof time for destructive superposition: Lb in (1) (dotted line), and for constructive
superposition: Lb in (2)(solid line).

2.5
2
1.5

One barTwo bar1

periodic, 3
0.5
0
0

20

40

60

80

Figure 2: Curves ampliﬁcation curve of surface elevation at the beach as functions of distace Lb for one-bar case (diamonds), two-bar (squares), periodic bar(triangles). All curves are periodically depending on Lb .

References[1] Pudjaprasetya, S.R. & Chendra, H.D., 2009 An Optimal Dimension of Submerged Parallel Bars as a Wave Reﬂector, Bulletin of Malaysian Mathematical
Soc., 2.
[2] Mei, C.C., Stiassnie, M., Dick K.-P. Yue, 2004, Theory and Applications of
Ocean Surface Waves, Advanced Series on Ocean Engineering Vol. 23, World
Scientiﬁc.
[3] Mattioli, F., 1990, Resonant Reﬂection of a Series of Submerged Breakwaters,
Il Nuovo Comento Vol. 13C, No.109
5, 823-833.
[4] Yu, J., Mei, C.C., 2000, Do longshore bars shelter the shore?, J. Fluid Mech.
404, 251-268,

Wave Energy Dissipation in Porous Media
Ikha Magdalena1),∗ , S.R. Pudjaprasetya2)
1, 2) Industrial & Financial Mathematics,
Fac. of Math. and Natural Sciences, Bandung Institute of Technology, Indonesia∗ e-mail: ikha.magdalena@yahoo.com

Abstract
An array of artificial reefs can be used as a submerged breakwater foroffshore protection. This type of breakwater provides environmental enhancement and aesthetics that are not found in conventional breakwaters.
An array of reefballs is assumed to be a submerged porous layer with certain
dimension and diffusive parameters. This research study the effectiveness ofa submerged porous media in dissipating wave energy of an incoming wave.
The approach is one-dimensional and linear.
The full governing equation for this problem consists of two Laplace
equations for a layer of ideal fluid, which is on top of a layer of fluid in
the submerged porous media. Boundary conditions are kinematic and dynamic along the surface, continuity of pressure and fluid velocity along theinterface, and impermeable bottom. Solving the full governing equation will
yield the following dispersion relation
εS tanh k(h1 − h2 ) cosh kh1 − αS sinh kh1
ω2
= S
gk
ε sinh kh1 tanh k(h1 − h2 ) − αS cosh kh1(1)

with αS = (C S − if S ). If we take parameter ω = 1, h1 = 1, h2 = 1.7,
ε = 0.4, C S = 0.4, f S = 1.5, dispersion relation (1) will results a wave
number k with negative imaginary part that results spatial dissipation ofwave amplitude, see Fig. 1 (Left). Fig. 1 (Middle & Right) show how|=(k)| depends on other variables.
From explicit formulas of velocity potential as solution of the full governing equation, we can derive an evolution equations in term of surface
elevation η and discharge Q:
1(nQ)x = 0n
Qt + c2 ηx = 0

ηt +

(2)
(3)

Coefficients n and c are complex numbers depending on variable bottom, andcharacteristic of permeable breakwater. We solve the equation numericallyusing finite volume method in a staggered grid. Numerical simulation shows
amplitude reduction of incoming waves, see Figure 2.

110

Figure 1: (Left) Solid line is the curve exp−i(kx−ωt) with k = 0.3117600956 −0.02450305720i, dotted line is its envelope. (Middle & Right) The curve of |=(k)|with respect to the non-dimensional variable d = (h2 − h1 )/h2 and r = ω 2 h2 /g.

Figure 2: (Left) Amplitude of incoming monochromatic wave reduces after passing
a submerged porous breakwater. (Right) Wave energy decreases due to submerged
porous breakwater.

References[1] Armono, H.D and Hall, K.R: Wave transmission on submerged breakwater made of hollow hemispherical shape artificial reefs, Proceeding of The
1st Coastal, Estuary and Offshore Engineering Speciality Conference, (2003),
CSC-185-1:10.
[2] Ching-Piao Tsai : Wave Transformation Over Submerged Permeable Breakwater on Porous Bottom, Ocean Engineering, (2006), 1623-1643.
[3] Dean, R.G., Dalrymple, R.A., Water Wave Mechanics for Engineers and Scientist, World Scientific, (1991).
[4] Stelling, G.S., Duinmeijer, S.P.A., A staggered conservative scheme for every
Froude number in rapidly varied shallow water flows, Int. J. Numer. Meth.
Fluids, (2003), 43, 1329-1354.

111

Multi-Compartment Modeling of Tumor Cells
Interacting
With Dynamic Chemotherapeutic Drug
Madhu Jain1)*, G. C. Sharma2)**
1) Indian Institute of Technology, Roorkee-247667(India)

2) Institute of Basic Science, Khandari, Agra-282002(India)*madhujain@sancharnet.in, **gokulchandra@sancharnet.in
Abstract

A tumor is an unusual growth of the tissues of the human body. Both benign
and malignant types of tumors require some kind of therapy to remove or reduce
their size. The growth of tumor is not regulated by normal body control
mechanisms. This unchecked growth is not beneficial to the organs of the body.
Chemotherapeutic drug has however been proved successful to great extent in the
reduction of the size of tumor. Benign tumor is a well-defined growth with smoothboundaries and is not dangerous in nature; it grows slowly in diameter andcompresses the nearby tissues. Malignant tumor usually has irregular boundaries
and attacks the surrounding area; it is poisonous in nature and takes the form ofcancer later on. Chemotherapy is a combination of two or more drugs whichadministrated to patients at different time and is very much useful in elimination or
reduction in the size of tumors.
The effect of drug on tumor growth can be mathematically modeled fromdifferent perspective. A three compartment mathematical model has beendeveloped to study the dynamics of tumor cells interacting with chemotherapeutic
drug. A system of differential equations that describe the dynamics of tumor cells,
is constructed. The model has been analyzed at equilibrium as well as transientstate. The steady state analysis has been done by using successive over relaxation
method whereas Runge Kutta method is employed for facilitating transition results.
The sensitivity analysis has also been performed to explore the effects of various
parameters on tumor cell population. The model analyzed demonstrate that the
drug with decaying factor has vital role in the decrease of proliferating cell
population in all proliferating compartment, however the drug has no effect on
non-proliferating cell population. It can be predicted through our model that
population of tumor cell can be controlled by adjusting the parameters during thecourse of chemotherapeutic drug administration. The outcomes of the investigation
may provide an insight for medical practitioners as well as pharmacists to getunderstanding the pharmacokinetic dynamics of the drug. Our model also reveals
that how equilibrium points and stability analysis can be employed in handling thetumor cell population.
112

On LaSalle’s Invariance Principle and Its
Application to Synchronize Hyperchaotic
Systems
Matthew Min-Hsiung Lin
Dept. of Mathematics, National Chung Cheng University∗ e-mail: mlin@math.ccu.edu.tw
Abstract
The notion of hyperchaotic synchronization is of fundamental importancebecause its wide range of important applications in secure communications,information processing, chemical reactions, and biological systems [1, 2, 3, 4].
In this talk, we want to focus on the adaptive control law to solve the synchronization problem. To synchronize two hyperchaotic systems, traditional
approaches often require the number of nonzero outputs of the controller
equal to the dimension of the state variables. However, in practice, it would
be desirable to reduce the number of outputs but also synchronize two hyperchaotic systems. Compared to traditional approaches, we show that that
two Rösssler hyperchaotic systems with known and unknown parameters canbe synchronized by only using three control functions. Theoretical proof isbased on the LaSalle’s invariant principle [5, 7], which is one of the very usefultheorems in dynamical systems and control theory. Plenty of applications arereferred to [6]. We believe that our idea of employing LaSalle’s invariant principle leads to a general strategy for claiming the synchronization between two
hyperchaotic systems. There is plenty of room for future research, includinga refinement of the control functions for more efficient synchronization anda generalization to more complex structures. We hope that our discussion
in this investigation offers a unified and an effectual avenue of synchronizing
between two given systems.

References[1] P. Badola, S. S. Tambe, and B. D. Kulkarni, Driving systems
with chaotic signals, Phys. Rev. A, 46 (1992), pp. 6735–6737.
[2] S. Bowong, Stability analysis for the synchronization of chaotic systems with different order: application to secure communications, Physics
Letters A, 326 (2004), pp. 102 – 113.

113

[3] G. Chen and X. Dong, From chaos to order: methodologies, perspectives, and applications, World scientific series on non-linear science:
Monographs and treatises, World Scientific, 1998.
[4] D. Ghosh, A. R. Chowdhury, and P. Saha, On the various kinds ofsynchronization in delayed duffing-van der pol system, Communicationsin Nonlinear Science and Numerical Simulation, 13 (2008), pp. 790 – 803.
[5] M. W. Hirsch, S. Smale, and R. Devaney, Differential Equations,Dynamical Systems, and an Introduction to Chaos - Second edition, 2ndedition, Academic Press, 2003.
[6] H. K. Khalil, Nonlinear Systems, 3rd edition, Prentice Hall, 2002.
[7] J. P. LaSalle, The Stability of Dynamical Systems, SIAM Press, 1976.

114

An unbiased estimator of noise correlations
under signal drift
Keiji Miura1,2),∗
1) PRESTO, Japan Science and Technology Agency
2) Department of Molecular and Cellular Biology, Harvard University∗ e-mail: kmiura@mcb.harvard.edu
Abstract

Cross correlations

Correlations in noises can play an important role in large-scale modelingand data analysis. In neuroscience, how reliably sensory information is represented by neural responses has been characterized by applying information
theoretic approach in a stochastic stimulus-response framework. It has beentheoretically shown that the correlations in response noises could be a major
determinant for coding capacities of sensory information by neurons [1] and,even in a very simple homogeneous network with tiny correlations, having
more neurons does not help at all [2]. Therefore, it is very important toestimate noise correlations accurately.
While signiﬁcant noise correlations has been reported
from almost all cortical areas,
it was pointed out that non0.1
stationarity such as a drift inmean responses could lead to
artiﬁcial correlations even if
0.05
there is no eﬀective correlation [3, 4, 5].
Although attempts to esti0
mate noise correlations underchanging environments have−20
020
been made, they worked only
Time shift
for speciﬁc cases [5].
For some statistical mod- Figure 1: Cross correlation function estiels, an unbiased estimator un- mated by the proposed method (full line) andder arbitrarily changing en- the conventional correlation coeﬃcient (dotvironments has been success- ted line) for the time shifted data. Σ11 =fully obtained in a simple, an- Σ22 = 1 and Σ12 = 0.1. Time dependent
alytically closed form by using signals ξ were generate by AR model.
the information geometry of
semiparametric models [6, 7].
115

In this paper I propose a semiparametric method for estimating noisecorrelations unbiasedly whatever the time course of mean drift is.
I consider a bivariate normal distribution for activities of two neurons:q(X; ξ, Σ−1 ) =1
1(X − ξ)′ Σ−1 (X − ξ)},exp{−
2π|Σ|1/2
2(1)where X = {x, y} and ξ = {ξ1 , ξ2 } are vectors. Here, I assume that the covariance matrix Σ is constant while the signals ξ can change over time. Especially, when the signals are randomly distributed, but, two consecutive signals
are the same, the distribution of activities {X} = {X1 , X2 } = {x1 , y1 , x2 , y2 }can be described as a mixture model:p({X}; Σ−1 , k) =∫k(ξ)q(X1 ; ξ, Σ−1 )q(X2 ; ξ, Σ−1 )dξ(2)

where k(ξ) denotes an unknown distribution of signals. I derive optimal−1−1−1
estimators of the three constant parameters Σ−1 = {Σ−1
11 , Σ12 (= Σ21 ), Σ22 }which work whatever k(ξ) is.

References[1] H. Sompolinsky, H. Yoon, K. Kang, M. Shamir: Population coding in neuronalsystems with correlated noise, Phys. Rev. E, 64 (2001), 051904.
[2] E. Zohary, M. N. Shadlen, W. T. Newsome: Correlated neuronal discharge
rate and its implications for psychophysical performance, Nature, 370 (1994),
140–35.
[3] A. S. Ecker, P. Berens, G. A. Keliris, M. Bethge, N. K. Logothetis, A. S.
Tolias: Decorrelated neuronal firing in cortical microcircuits, Science, 327(2010), 584–7.
[4] A. Renart, J. de la Rocha, P. Bartho, L. Hollender, N. Parga, A. Reyes, K.
D. Harris: The asynchronous state in cortical circuits, Science, 327 (2010),
587–90.
[5] W. Bair, E. Zohary, W. T. Newsome: Correlated firing in macaque visualarea MT: time scales and relationship to behavior, J Neurosci., 21 (2001),
1676–97.
[6] S. Amari, M. Kawanabe: Information geometry of estimating functions insemi-parametric statistical models, Bernoulli, 3 (1997), 29–54.
[7] K. Miura, M. Okada, S. Amari: Estimating spiking irregularities under changing environments, Neural Comput., 18 (2006), 2359–86.

116

Construction of Probabilistic Boolean
Networks: A Maximum Entropy Rate
Approach
Wai-Ki Ching1),∗ ,1) Advanced Modeling and Applied Computing Laboratory
Department of Mathematics
The University of Hong Kong
Pokfulam Road, Hong Kong.
∗ e-mail: wching@hku.hk
Abstract
Modeling genetic regulatory networks is an important problem in genomic research. Boolean Networks (BNs) and their extensions Probabilistic
Boolean Networks (PBNs) have been proposed for modeling genetic regulatory interactions.
In a PBN, its steady-state distribution gives very important information
about the long-run behavior of the whole network. However, one is alsointerested in system synthesis which requires the construction of networks.
The inverse problem is ill-posed and challenging, because there may be many
networks or no network having the given properties, and the size of the problem is huge. The construction of PBNs from a given transition-probability
matrix and a given set of BNs is an inverse problem of huge size.
Here we propose a maximum entropy approach for the above problem.
Newton’s method in conjunction with the Conjugate Gradient (CG) method
is then applied to solving the inverse problem. We investigate the convergence rate of the proposed method. Numerical examples are also given to
demonstrate the eﬀectiveness of our proposed method.
This is a joint work with X. Chen, X.S. Chen, Y. Cong and N. Tsing [3].

References[1] W. Ching, X. Chen and N. Tsing, Generating Probabilistic Boolean Networks from a Prescribed Transition Probability Matrix, IET on Systems
Biology, 6 (2009) 453-464.
[2] S. Zhang, W. Ching, X. Chen and N. Tsing, Generating Probabilistic
Boolean Networks from a Prescribed Stationary Distribution, Information Sciences, 180 (2010) 25602570.
117

[3] X. Chen, W. Ching, X.S. Chen, Y. Cong and N. Tsing, Construction of
Probabilistic Boolean Networks from a Prescribed Transition Probability
Matrix: A Maximum Entropy Rate Approach, to appear in East Asian
Journal of Applied Mathematics, (2011).

118

Compensation Technique for Transmission Delay andPacket Loss in Networked Control System based on
Lagrange Interpolation1),*

Inseok Yang

2),

2)

, Donggil Kim , Dongik Lee1) Department of Industrial Applied Mathematics, Kyungpook National University,
Korea
2) School of Electronics Engineering, Kyungpook National University, Korea* e-mail: hllover@knu.ac.kr
Abstract
In this paper, the compensation technique for transmission delay/packet loss in
networked control system (NCS) is proposed. NCS is a feedback control system
wherein the control loops are closed through a communication network thanks torecent advances in microelectronics [1]. NCS has many advantages compared to
traditional control systems, including ease to maintenance, reduced system space
and weight by replacing mechanical linkages to network nodes, etc.
However, by inserting communication network in feedback control loop, newproblems including transmission delay/packet loss among network nodes occur asshown in Fig. 1. These problems degrade the system performance and even makesystems unstable. To overcome these problems, this paper proposes a compensation
technique based on Lagrange polynomial interpolation [2]. The proposed Lagrange
polynomial compensation method can compensate the lost data by updating thepast received data. For this reason, the compensation module can be implemented
as an independent network device, as shown in Fig. 2. The equation of the proposed
Lagrange polynomial compensation method can be expressed as follow: fˆ (nT ) =2

 f ((n − 3)T )  3

nT − (n − j )T f ((n − i)T ) = [1, − 3, 3]  f ((n − 2)T )  aT vn (1)
∑
∏


(
n
i
)
T
(
n
j
)
T
−
−
−
i =1
j =1
 j ≠i f ((n − 1)T ) 3

where, T is a sampling period. In (1), the lost data can be compensated by updating
vn in real time.
The performance of the proposed technique is evaluated by numericalsimulation results applied to an “innovative control effectors (ICE)” aircraft [3].
The aircraft maneuvers with maximally 3.5 deg of error between the normal modeland the uncompensated model as shown in Fig. 3 and 4. However, if the aircraft
adopts the proposed compensation module, then the error is reduced significantly,
to less than 0.5 deg. Hence, the simulation results show that the proposed method
can compensate lost data fast and effectively.

119

References[1] Y. Tipsuwan, M.-Y. Chow: Control Methodologies in Networked Control
Systems, Control Engineering Practice, 11 (2003), 1099--1111.
[2] S. C. Chapra, R. P. Canale: Numerical Methods for Engineers, McGrawHill,(2006).
[3] S. R. Wells: Application of Sliding Mode Methods to the Design ofReconfigurable Flight Control Systems, Ph. D. Dissertation, University ofCalifornia, (2002).
System coordinatey (t )

Plant

Actuator

TDC

Sensory (kT )

u (kT − τ ca )

τ caNetworkτ cay (kT − τ sControlleru (kT )

Actuatory (kT )
τ sc

Networky (kT − τ sc )Control coordinate

Controller

Fig. 1. The structure of networked control
systems.

Sensoruˆ(kT )

u (kT − τ ca )

τ scu (kT )y (t )Plantyˆ( kT )TDC

Fig. 2. The proposed architecture of networked
control system.

Fig. 3. The results of turn reversal maneuver of the ICE aircraft with and without compensation.

Fig. 4. The errors of turn reversal maneuver of the ICE aircraft with and without compensation.

120

Mathematics of Crime
Andrea BERTOZZIUniversity of California Los Angeles* e-mail: bertozzi@math.ucla.edu

Abstract
There is an extensive applied mathematics literature developed for problems in the
biological and physical sciences. Our understanding of social science problems
from a mathematical standpoint is less developed, but also presents some very
interesting problems, especially for young researchers. This lecture uses crime as acase study for using applied mathematical techniques in a social science applicationand covers a variety of mathematical methods that are applicable to such problems.
We will review recent work on agent based models, methods in linear and nonlinear
partial differential equations, variational methods for inverse problems andstatistical point process models. From an application standpoint we will look at
problems in residential burglaries and gang crimes. Examples will consider both``bottom up'' and ``top down'' approaches to understanding the mathematics ofcrime, and how the two approaches could converge to a unifying theory.

121

Numerical Studies of Lorentz Transformation*

David C. Ni

Direxion Technology
9F, No. 177-1, Ho-Ping East Rd., Section 1, Taipei, Taiwan* e-mail: davidcni@yahoo.com
Abstract
The invariance of Lorentz transformation, which the relativity theories are based
upon, has been under extensive tests and elaborated extensions in several areas
recently. In the area of astrophysics, the origin of magnetic fields in conjunction
with generalized vorticity is connected to the space-time distortion [4]. In the areaof quantum gravity, the speed of light is predicted to vary at Planck scale [1 andreferences therein]. At the high-energy or ultraviolet limit, new developments, suchas lattice gauge theory, which evolving from Yang-Mills gauge theory [8] andrenormalization [2, 3, 7], still have debated issues as those in quantum gravity. We
present numerical studies of Lorentz transformation with focus on the nonlinear
extensions of the equations. We find interesting results such as fractal patterns,
chaos, nonlinear-to-linear degeneracy, continuum-to-discreteness transition,
vorticity and symmetry broken.
The efforts in testing Lorentz invariance are proposed with the new theoreticalmodels based on phenomenological modifications of the usual Lorentz invariant
dispersion law. One of the models modifies the dispersion law E2 = p2c2 + m2c4 byadding nonlinear terms in powers of particle momentum as follows:
E2 = p2c2 + m2c4 + h1p + h2p2 + h3p3 + …

(1)

Where c is the light velocity, E and p are the energy and momentum respectively, mis the mass of particle, and hj is the coefficient. Since these proposals are still
conjectural, we are motivated to examine the dispersion law and related
phenomenological modifications by numerical analysis [5, 6].
Given two inertial frames with different velocities, u and v, the observed velocity,u’ , from v-frame is as follows:2u’ = ( u - v ) / ( 1 - vu/c )(2)

We set c2 = 1 and the multiply a phase connection, exp(iχ(u)), to the normalized
complex form of the equation (2) based on the concept of gauge transformation [3]as follows:122f = (u’/u) = exp(iχ(u))(1/u)[(u-v)/( 1-uv)](3)

We adopt the appraoch by iterating the function f for all the points on the complex
plane and examine the convergent subsets. We find several interesting observations
as follows:(1) The convergent subsets show fractal patterns(2) For particular v values, chaotic phenomena occurs(3) Nonlinear-to-linear degeneracy occurs when v value, which corresponding to
the normalized velocity of observing frame, is approaching to unity.
(4) Abrupt continuum-to-discreteness transition when v value further approaching
to unity in the range of 10-16.
(5) After the transition occurs, the convergent subsets are symmetrical to both realand imaginary axes.
References[1] A. A. Abdo et al, A limit on the variation of the speed of light arising fromquantum gravity effects, Nature 462 (Nov., 2009) 331-334[2] G. 't Hooft and M. Veltman, Regulation and renormalization of gauge fields,
Nuclear Physics B 44, (1972) 189–219[3] G. 't Hooft, ed., 50 Years of Yang-Mills theory, (2005) World Scientific[4] S. M. Mahajan and Z. Yoshida, Twisting Space-Time: Relativistic Origin of
Seed Magnetic Field and Vorticity, Phys. Rev. Lett. 105 (2010) 095005[5] David C. Ni and C. H. Chin, Classification on Herman Rings of Extended
Blaschke Equations, Differential Equations and Control processes, #2,
Article1(2010)( http://www.neva.ru/journal/j/EN/numbers/2010.2/issue.html)[6] D. C. Ni, N. B. Ampilova, B. Alina, and C. H. Chin, Chaotic Behavior andBifurcation in N-Body Systems, Proceedings for 8th AIMS International
Conference on Dynamical Systems, Differential Equations, and Applications(2011)[7] Kenneth G. Wilson, The renormalization group: critical phenomena and the
Kondo problem, Rev. Mod. Phys., 47 (1975), 773-839[8] C. N. Yang and R. Mills, Conservation of Isotopic Spin and Isotopic Gauge
Invariance, Physical Review 96 (1954) 191–195

123

Numerical studies on the Monge-Ampéreequation arising from free-form design ofgeometric optics
Yu-Lin Tsai1) , Ming-Chiang Jiang2) and Chin-Tien Wu2),∗
1) Department of Photonics and Institute of Electro-Optical Engineering,
Chiao-Tung University, Taiwan.
2) Department of Applied Mathematics, Chiao-Tung University, Taiwan.
∗ e-mail: ctw@math.nctu.edu.tw
Abstract
The Monge-Ampere (MA) equation is a fully nonlinear PDE arising invarious research fields including differential geometry, mass transportation,
geostrophic fluid and optical design. Numerical methods such as finite difference methods proposed by Oberman etc [7], Dean and Glowinshi [3, 4] andfinite element methods proposed by Feng and Neilan [5, 6] and Awanou [8]have been successfully applied to solve the standard MA equation. In thistalk, we introduce some of known theoretical results obtained by Caffarelli[9], Oliker [1], Mader [12] and Wang [11] etc., and the above mentioned numerical methods. Particularly, we employee the vanished-moment method tosolve the MA equation arising from the free form surface (FFS) design problem in geometric optics[2] [13]. To ensure the FFS obtained from numericalcomputation is continuous and to be able to compute the curvature of the
FFS easily, we solve the linearized MA problem using BCIZ finite element
method[10]. The discrete linear system is solved by the conjugate gradientmethod (CG) in which multigrid approximation of the discrete biharmonic
operator is used as the preconditioner of the CG method. Accuracy and robustness of our approach are demonstrated on several benchmark examples.
We compare our numerical results with the results obtained by Feng and
Neilan, Froese and Oberman, and Dean and Glowinski.

References[1] V. Oliker, Mathematical aspects of design of beam shaping surface ingeometrical optics., In: Trends in Nonlinear Analysis, Springer-Verlag,
pp. 191–222, 2002.
[2] J. S. Schruben, Formulation of a reflector-design problem for a lighting fixture., J. Opt. Soc. Am. 62, pp. 1498– 1501, 1972.

124

[3] E. J. Dean and R. Glowinshi, Numerical solution of the twodimensional elliptic Monge-Ampére equation with Dirichlet boundary
conditions: an augmented Lagrangian approach., C. R. Acad. Sci. Paris,
Ser. I 336, pp. 779–784, 2003.
[4] E. J. Dean and R. Glowinshi, Numerical methods for fully nonlinear
elliptic equations of the Monge-Ampére type., Comput. Methods Appl.
Mech. Engrg. 195, pp. 1344–1386, 2006.
[5] X. Feng and M. Neilan, Mixed finite element methods for the
fully nonlinear Monge-Ampere equation based on the vanishing moment
method., SIAM J. Numer. Anal., 47, pp. 1226-1250, 2009.
[6] X. Feng and M. Neilan, Vanishing moment method and moment
solutions for fully nonlinear second order partial differential equations.,
J. Sci. Comput., 38, pp. 74-98, 2009.
[7] J. D. Benamou, B. D. Froese and A. D. Oberman, Two numerical methods for the elliptic Monge-Ampere equation., ESAIM: Mathe.
Modell. Numer. Anal., pp. 737–758, 2009.
[8] G. Awanou, Numerical methods for fully nonlinear elliptic equations.,
Inter. Conf. App. Thm., 2010.
[9] L. A. Caffarelli and X. Cabre, Fully nonlinear elliptic equations.,
American Mathematical Society Colloquium Publications. 43, American
Mathematical Society, Providence, RI, 1995.
[10] G. P. Bazeley, Y. K. Cheung, B. M. Irons and O. C.
Zienkiewicz, Triangular elements in plate bendingconforming and nonconforming solutions, Proc. Conf. on Matrix Methods in Structural Mechanics, WPAFB, Ohio, pp. 547-576, 1965.
[11] P.Guan and X. J. Wang, On a Monge-Ampere equation arising ingeometric optics J. Differential Geometry, 48, pp. 205-223, 1998.
[12] L. Marder,Uniquessness in reflector mappings and the Monge-Ampere
equation Proceeding of Royal Society of London. Series A, Mathematicaland Physical Sciences, 378, No. 1775, pp.529-538, 1981.
[13] P. Benitez, J. C. Minano, J. Blen, R. Mohedano, J. Chaves,
O. Dross, M. Hernandez andW. Falicoff Simultaneous multiple
surface optical design methof in three dimensions., Opt. Eng. 43(7). pp.
1489-1502, 2004.

125

Renormalization approach to solve
Langevin’s equations
Takashi Ichinomiya
School of Medicine, Gifu University
e-mail: tk1miya@gifu-u.ac.jp
Abstract
The importance of molecular dynamics(MD) simulations in physics, chemistry and biology has been increased in recent years. Unfortunately, when wewant to investigate very rare events, such as chemical reactions, the computational costs of the MD can be extremely large. Several methods have been
proposed to overcome this diﬃculty[1, 2, 3].
In this talk, we propose a new method to accelerate MD simulation, using the idea of renormalization. The application of renormalization to MD
simulation has been proposed by Faccioli[4]. However our method use the different renormalization technique, inspired by the work of Chen, Goldenfeld,
and Oono[5].
In the following, we consider the following Langevin’s equation√dxi= f (x) + Dξi (t),dt(1)

where x = (x1 , x2 , · · · , xn ), and ξi (t) is Gaussian white noise, which satisﬁeshξi (t)ξj (t0 )i = δ(t − t0 )δij . If we use the Euler-Maruyama scheme to solve thisequation numerically, x(∆t) is given by√xi (∆t) = ∆tf (x(0)) + D∆tηi (t),(2)
where ηi is Gaussian random noise. Therefore the distribution of x(∆t) is
given asP (x(∆t)|x0 ) ∝ exp[−(x − ∆tf (x0 ))T A0 (x − ∆tf (x0 ))],

(3)

where (A0 )ij = δij (2D∆t)−1 .
The Euler-Maruyama scheme is valid only when ∆t is small enough. Ourbasic idea is to assume that P (x(t)|x0 ) can be approximatedP (x(t)|x0 ) ∝ exp[−(x − tf (x̃(t)))T At (x − tf (x̃(t)))],

126

(4)

for larger t, where x̃(t) and At is renormalization parameters. To calculatex̃(t) and At , we use the equation∫

P (x(2t)|x0 ) =dx0 P (x0 (t)|x0 )P (x(t)|x0 ).

(5)

By the Taylor expansion of f (x) around x = x̃(t), we get the approximatedequation for x̃(2t) and A2t asx̃(2t) = x̃(t) + t2 f1 (x̃(t))f (x̃(t)) + (x̃(t) − x0 )
and
A2t =

A t
− (At f1 (x̃(t)) + f1 (x̃(t))At ),
2
4

(6)

(7)

∂f
where f1 = ∂x
. By applying these equations n-times,we can obtain P (x(2n ∆t)|x0 ).
Eqs.(6) and (7) serves us a new simulation algorithm. In my talk, we will
present our new method in detail, and discuss on the advantages and disadvantages of this method, with some results of numerical simulations.

References[1] C. D. Dellago, P. G. Bolhuls, F. S. Csajka and D. Chandler, Transition
path sampling and the calculation of rate constants, The Journal ofChemical Physics, 108(1998) 1964-1977.
[2] A. F. Voter, A method for accelerating the molecular dynamics simulation of infrequent events, Journal of Chemical Physics, 106(1997),
4665-4677.
[3] G. Hummer and I. G. Keverkidis, Coarse molecular dynamics of a peptide fragment: Free energy, kinetics, and long-time dynamics computations, Journal of Chemical Physic, 118(2003) 10762-10773.
[4] P. Faccioli, Molecular dynamics at low time resolution, The Journal ofChemical Physics, 133(2010), 164106.
[5] L-Y. Chen, N. Goldenfeld, Y. Oono, Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive
perturbation theory, Physical Review E, 54(1996) 376-394.

127

Consistent approximation-sampling withmulti pre and post filtering
Kil Hyun Kwon1),∗ , Jaekyu Lee1)
1) Department of Mathematical Sciences, KAIST, S. Korea.
∗ e-mail: khkwon@kaist.edu
Abstract
Sampling is the process of representing signals of a continuous variable t
by sequence of numbers, that is, discrete signal representation.
If a signal f (t) is band-limited, that is, f (t) ∈ P Wπ∑, then f (t) is completely characterized by its samples {f (n)}n∈Z as f (t) =f (n)sinc(t−n) =n∈Z(c ∗ φ)(t), where c = {f (n)}n∈Z ∈ l2 and φ(t) = sinc(t) is the interpolation
kernel(or interpolator).
If a signal f (t) ∈ L2 (R) is not band-limited, then the standard procedure
is to apply an additional prefiltering with the ideal low pass filter χ[−π,π] (ξ) =bφ(ξ),φ(t) = sinc(t) before sampling. Then the corresponding band-limited∑
reconstruction fe(t) =hf (t), φ(t − k)iφ(t − k) is the orthogonal projectionk∈Zof f (t) onto P Wπ .
More generally for any Riesz generator φ(t) ∈ L2 (R),∑e − k)iφ(t − k)hf (t), φ(tfe(t) =

(0.1)

k∈Z

e is theis the orthogonal projection of f (t) ∈ L2 (R) onto V (φ), where φ(t)dual Riesz generator of φ(t).
Note that fe(t) in (0.1) is the L2 -least square approximation of f ∈ L2 (R)ein V (φ). We call h(t) := φ(−t)
the analysis function and f 7−→ fe as in (0.1)
the ”approximation-sampling” [1].
However, in most applications, the analysis function is given a priori as
the impulse response h(t) of the acquisition device. So, [6] introduces theconsistent approximation-sampling procedure, where the consistency meansthat the original signal and its approximation are essentially equivalent toobservers(through sensors or acquisition devices). Many authors[2, 3, 4, 5,
7, 8] have extended it several ways.
In this work, we consider the consistent approximation-sampling procedure satisfying m-shift-invariance when the samples hf (t), ψi (t−qk)i (1 ≤ i ≤M, k ∈ Z) of f (t) ∈ L2 (R) and the reconstruction filters {φj (t)|1 ≤ j ≤ N }128

are given, where {ψi (t − qk)|1 ≤ i ≤ M, k ∈ Z} is a Riesz sequence ,{φj (t)|1 ≤ j ≤ N } is a Riesz generator and q is a rational number. Then wefind some equivalent conditions that the approximation exists. Also, we discuss the properties of the approximation and which approximation is closest
to the least square solution (0.1).

References[1] A. Aldroubi and M. Unser: Sampling procedure in function spaces andasymptotic equivalence with Shannon’s sampling theory, Numer. Func.
Anal. Opt., 15 (1994), 1–21.
[2] Y. Eldar: Sampling with arbitary sampling and reconstruction spaces
and oblique dual frame vectors, J. Fourier Anal. Appl., 9 (2003), 77–96.
[3] Y. Eldar and T. G. Dvorkind: A minimum squared-error framework forgeneralized sampling, IEEE T. Signal Proces., 54 (2006), 2155–2167.
[4] Y. Eldar and T. Werther: Generalized framework for consistent sampling in Hilbert spaces, Int. J. Wavelets. Multi., 3 (2005), 347–359.
[5] A. Hirabayashi and M. Unser: Consistent sampling and signal recovery,
IEEE T. Signal Proces., 55 (2007), 4104–4115.
[6] M. Unser and A. Aldroubi: A general sampling theory for nonideal
acquisition devices, IEEE T. Signal Proces., 42 (1994), 2915–2925.
[7] M. Unser and J. Zerubia: Generalized sampling: Stability and performance analysis, IEEE T. Signal Proces., 45 (1997), 2941–2950.
[8] M. Unser, J. Zerubia: A generalized sampling theory without bandlimiting constraints, IEEE T. Circuits-II., 45 (1998), 959–969.

129

A Tailored Finite Point Method forConvection Diffusion Reaction Problems with
Variable Coefficients
Yintzer Shih1,∗ , R. Bruce Kellogg2
1. Department of Applied Mathematics, National Chung Hsing University
Taichung 402, Taiwan
2. Department of Mathematics, University of South Carolina, SC 29208,
USA∗ e-mail: yintzer shih@email.nchu.edu.tw
Abstract
We present a tailored finite point method (TFPM) [1, 2, 3] for solving theconvection-diffusion-reaction equation with variable coefficients. We develop
an algorithm to construct a streamline-aligned grid for the TFPM. Our numerical tests show for small diffusion coefficient the TFPM solution resolves
the internal and boundary layers regardless the mesh size, and depicts that
TFPM method with a streamline grid has excellent performance comparedwith the tailored finite point method and a streamline upwind finite element
method when ε is small.

References[1] H. Han, Z. Huang and R.B. Kellogg: The tailored finite point method
and a problem of P. Hemker, Proceedings of the International Conference on Boundary and Interior Layers - Computational and Asymptotic
Methods, Limerick, July 2008.
[2] Y. Shih, R.B. Kellogg and P. Tsai: A Tailored Finite Point Method forConvection-Diffusion-Reaction Problems, J. of Sci. Compu., 43 (2010),
pp. 239-260.
[3] Y.T. Shih, R.B. Kellogg and Y. Chang: Characteristic Tailored Finite
Point Method for Convection-Dominated Convection-Diffusion-Reaction
Problems, J. of Sci. Compu., DOI: 10.1007/s10915-010-9433-7.

130

On the interpolation constants overtriangular elements
Kenta Kobayashi
Institute of Science and Engineering, Kanazawa University
e-mail: Kenta.k@staﬀ.kanazawa-u.ac.jp
Abstract
We present new formulas which give sharp upper bounds for the interpolation constants on the triangles. Our formulas provide better upper bounds
than the former ones. Furthermore, from the point of practical use, our formulas are convenient and easy to calculate. In the proof of the formulas, weemploy numerical veriﬁcation method.
Let T be the triangle in R2 and V 1,1 (T ), V 1,2 (T ), V 2 (T ) be the function
spaces deﬁned by
½
¾
¯ Z
¯
1,1
1
V (T ) = φ ∈ H (T ) ¯
φ dxdy = 0 ,
T
½
¾
¯ Z
¯
1,2
1
V (T ) = φ ∈ H (T ) ¯
φ ds = 0, k = 1, 2, 3 ,
γk
¯
n
o
¯
V 2 (T ) = φ ∈ H 2 (T ) ¯ φ(pk ) = 0, k = 1, 2, 3 ,where p1 , p2 , p3 and γ1 , γ2 , γ3 are the vertices and sides of T respectively.
Then, it is known that the following constants C1 (T ), C2 (T ), C3 (T ), C4 (T )exist:
C1 (T ) =∥φ∥L2 (T ),φ∈V 1,1 (T )\0 ∥∇φ∥L2 (T )C2 (T ) =∥φ∥L2 (T ),φ∈V 1,2 (T )\0 ∥∇φ∥L2 (T )C3 (T ) =∥φ∥L2 (T ),φ∈V 2 (T )\0 |φ|H 2 (T )C4 (T ) =∥∇φ∥L2 (T ),φ∈V 2 (T )\0 |φ|H 2 (T )

sup

sup

sup

sup

where |φ|H 2 (T ) is a H 2 semi-norm of φ deﬁned by|φ|2H 2 (T ) = ∥φxx ∥2L2 (T ) + 2∥φxy ∥2L2 (T ) + ∥φyy ∥2L2 (T ) .
For these constants, we have obtained the formulas which give sharp upperbound of Cj (T ) asCj (T ) < Kj (T ),j = 1, 2, 3, 4.

Concrete form of the formulas Kj (T ) will be shown in the talk.
131

A posteriori estimates of inverse linear
ordinary diﬀerential operators
T. Kinoshita1),∗ , T. Kimura2) and M.T. Nakao2)
1) Research Institute for Mathematical Sciences, Kyoto University
2) Sasebo National College of Technology∗ kinosita@kurims.kyoto-u.ac.jp
Abstract
We present constructive a posteriori estimates of inverse operators for
initial value problems in linear ordinary diﬀerential equations on a boundedinterval. Here, “constructive” indicates that we can obtain the bounds of theoperator norm in which all constants are explicitly given or are represented
in a numerically computable form. In general, it is diﬃcult for us to obtain
this estimates of this inverse operators in the meaning of a priori. We propose the technique for obtaining the a posteriori estimates by using Galerkinapproximation of inverse operators.

References[1] M.T. Nakao, N. Yamamoto and S. Kimura, On the Best Constant in theError Bound for the H01 -Projection into Piecewise Polynomial Spaces,
Journal of Approximation Theory, 93 (1998), 491–500.
[2] M.T. Nakao, K. Hashimoto, Y. Watanabe, A numerical method to verify
the invertibility of linear elliptic operators with applications to nonlinear
problems, Computing, 75 (2005), 1–14.
[3] M.T. Nakao and K. Hashimoto, Guaranteed error bounds for ﬁnite element approximations of noncoercive elliptic problems and their applications, Journal of Computational and Applied Mathematics, 218 (2008),
no. 1, 106–115.

132

Numerical veriﬁcation for solution existence ofsemilinear elliptic partial diﬀerential equation on
arbitrary polygonal domain
Akitoshi Takayasu 1) , Xuefeng Liu 2),∗ , Shin’ichi Oishi 2)
1) Graduate School of Fundamental Science and Engineering, Waseda University
2) Faculty of Science & Engineering, Waseda University; CREST/JST∗ e-mail: xﬂiu@aoni.waseda.jp
Abstract
We consider the semilinear elliptic partial diﬀerential equation on polygonal domain of arbitrary shape.
−∆u = f (u) on Ω, associated with boundary condition .

(1)

To verify the existence of solutions to equation (1), there have appeared
several methods adopting invariant point theories, c.f. [2, 3, 4], which declarethe solution existence in a small range nearby approximate solution.
In case that the domain has non-convex corner, the solution u of equation(1) can be singular one, that is, the L2 integration of second order derivativeof u is not bounded, which brings diﬃculty for constructing error estimation.
By using the Prager-Synge’s theorem, we successfully developed computablea priori error estimation as below [1]:k∇(u − uh )kL2 (Ω) ≤ Ch kf k2 ,where uh is approximate solution from linear conforming FEM method and
Ch is computable quantity depending only on domain and mesh size.
By further applying the Newton-Kantorovich theorem, see also [4], we
can verify the solution existence of equation (1) even on polygonal domainof arbitrary shape.
In the presentation, we will display numerical computation results todemonstrate the eﬀciency of proposed method.

References[1] Xuefeng Liu and Shin’ichi Oishi, Veriﬁed eigenvalue evaluation for thelaplacian on arbitrary polygonal domain, submitted to Numerische
Mathematik, Dec,2010[2] Mitsuhiro T. Nakao and Yoshitaka Watanabe, Numerical veriﬁcation
methods for solutions of semilinear elliptic boundary value problems,Nonlinear Theory and Its Applications, IEICE, vol.2, no.1, pp.2-31.,
2011

133

[3] Plum.M, Computer-assisted enclosure methods for elliptic diﬀerential
equations, Linear Algebra and its Applications, Vol.324, pp.147-187.,
2001[4] Shin’ichi Oishi, Akitoshi Takayasu and Takayuki Kubo, Numerical veriﬁcation of existence for solutions to Dirichelt boundary value problemsof semilinear elliptic equations, Preprint submited to Journal of Computational & Applied Mathematics.

134

Reconstruction of the density oftree-shaped networks from boundary
measurements of waves
WADA Naoki1),∗
1) Graduate School of Informatics, Kyoto University∗ e-mail: naoki@acs.i.kyoto-u.ac.jp
Abstract
On finite graphs whose edges are identified with open intervals, we consider the inverse problem for wave equation to determine the unknown coeﬃcient from boundary measurements. In general case this inverse problem
is not uniquely solvable. We restrict ourselves in the case of trees and givea reconstruction by reducing the inverse problem to an infinite family ofFredholm integral equations of the second kind.

1

Waves on Graphs

Let Ω be a finite connected tree, which is consisting of edges E = {ej }K
j=1
and vertices V = {vj }Lj=1 . Each edge ej is identified with an open interval
and parameterized through this identification. The vertices with more than
one incident edge are called the interior vertices and the others are called theboundary vertices. The set of the interior vertices and the boundary vertices
are denoted by N and Γ, respectively.
Let ρ( · ) be a positive C 2 -smooth function on Ω\V and let T be a positive
number. We consider the wave equation
ρ(x)utt (t, x) − uxx (t, x) = 0, t ∈ (0, T ), x ∈ Ω \ Vu(t, · ) ∈ C(Ω), t ∈ [0, T ]∑∂ν u|ej (t, v) = 0, t ∈ [0, T ], v ∈ N

(1)
(2)
(3)

ej ≻v

u(0, x) = ut (0, x) = 0, x ∈ Ωu(t, γ) = f (t, γ), t ∈ [0, T ], γ ∈ Γ,

(4)
(5)

where u|ej is the restriction of u on the edge ei , ∂ν is the outward normal at
endpoints of intervals and the sum in (3) is taken all over the edges incident
to the interior vertex v ∈ N , which condition is called Kirchhoﬀ condition.
135

For this initial boundary value problem, we can difine the operator RT whichmaps Dirichlet data uf |[0,T ]×Γ = f to(the Neumanndata( ∂ν uf |[0,T ]×Γ).
)
We difine the operator C T : L2 (0, T ) × Γ → L2 (0, T ) × Γ by the
relation∫( T
C f, g)L2 ((0,T )×Γ) =uf (T, x) ug (T, x) ρ(x) dx.
Ω

This operator C T is expressed through the response operator RT [2].

2Reconstruction

We give a local reconstruction of ρ from the boundary using the boundary
measurements RT . On an edge e which connected to a boundary vertex
γ0 , we introduce some formulae [1]. Let α, β T be functions whose supportsT
are contained in [0, T ] × {γ0 } and α(t,( Tγ0 )) = t, β (t, γ0 ) = T − t. For thesefunctions, we have the relation lim R α (t) = −ρ(γ) and equations
t→+0

C T uf = (RT )∗ β T(6)

C T ug = β T .

(7)Functions ft0 , gt0 that satisfies supp ft0 , supp gt0 ⊂ [T − t0 , T ] × {γ0 } and (6),(7) in [T − t0 , T ] × {γ0 } exist uniquely respectively. The equations (6),(7)are expressed throgh the boundary measurements RT and have the formof Fredholm integral equations of the second kind. These solution satisfies
) 41
)1
(
(
ρ(x̂(t0 )) 4
0 ))
,
lim
g
x̂(t0 ), where(t)=that limt→t0 +0 ft0 (t) = ρ(x̂(tt→t0 +0 t0
ρ(γ0 )
ρ(γ0 )∫x √x̂(t) is the inverse function of τ (x) = γ0 ρ(x) dx. τ (x) is the so-called
optical distance, which is the travel time for waves from γ0 to x. Using these
relations, by solving the equations (6)(7), we can reconstruct the unknown
coeﬃcients ρ on the edge e.
By generalizing these process mentioned above, we show a reconstructionof the unknown coeﬃcients on trees.

References[1] M. Belishev: Boundary control and inverse problems: The onedimansional variant of the BC-method, J. Math. Science, Vol.155, No.3,(2008) pp.343–378[2] M. Belishev, N. Wada: On revealing graph cycles via boundary measurements, Inverse Problems, 25 (2009)136

Cubic helical splines with Frenet-frame
continuity
Chang Yong Han1),∗ Song-Hwa Kwon2)
1) Department of Applied Mathematics, Kyung Hee University2) Department of Mathematics, The Catholic University of Korea∗ e-mail: cyhan@ymail.com
Abstract
A parametric curve r(t) in R3 is called a helix if its tangent vector maintains a constant angle against a fixed axis [7]. A helix admits polynomial
parametrization only if it satisfies Pythagorean-hodograph (PH) condition[4], i.e., its speed is a polynomial in the curve parameter [1]. Although helical PH curves comprise a proper subset of spatial PH curves, all spatial PH
cubics are helical [2, 3]. As the simplest class of PH curves and polynomial
helices, spatial PH cubics have been extensively studied in the context ofHermite interpolation [9, 5, 8, 6].
Our focus in this presentation is on the geometric spline interpolating asequence of points p0 , . . . , pn in R3 . This subject has also been discussed
in the aforementioned papers as a straightforward extension of G1 Hermite
interpolation problem: at each pk a unit tangent tk is either derived from theoriginal curve [5], or assigned according to a shape-preserving interpolation
scheme [8] to generate a sequence of PH cubic segments r1 (t), . . . , rn (t) that
constitute a G1 spline for the given data points. Instead of directly specifying
the unit tangents, we determine each tk in such a way that the Frenet framesof the adjacent PH cubic segments are coincident at their juncture point
pk , therefore improving the continuity of the resulting spline to be not only
tangent-continuous but also Frenet-frame continuous. It turns out that toensure a finite number of solutions, an additional condition is needed foreach segment. We require each rk (t) to satisfy a certain symmetry property
which is geometrically attractive and facilitates the subsequent analysis. The
construction of rk (t) involves a new kind of geometric Hermite data consistingof an initial Frenet frame and the displacement vector ∆pk = pk −pk−1 . The
existence condition of a PH cubic interpolant is formulated in terms of theposition of ∆pk relative to the Frenet frame. Once an initial Frenet frame at
p0 is provided, the entire collection of PH cubic segments can be sequentiallygenerated using the Frenet frame of rk−1 (t) at pk−1 as the initial Hermite
data for rk (t).

137

References[1] R.T. Farouki: Pythagorean-Hodograph Curves: Algebra and Geometry
Inseparable, Springer, Berlin, 2008.
[2] R.T. Farouki, C. Giannelli, A. Sestini: Helical polynomial curves anddouble Pythagorean hodographs I. Quaternion and Hopf map representations. Journal of Symbolic Computation, 44 (2009), 161–179.
[3] R.T. Farouki, C. Giannelli, A. Sestini: Helical polynomial curves anddouble Pythagorean hodographs II. Enumeration of low-degree curves.
Journal of Symbolic Computation, 44 (2009), 307–332.
[4] R.T. Farouki, C.Y. Han, C. Manni, A. Sestini: Characterization andconstruction of helical polynomial space curves. Journal of Computational and Applied Mathematics, 162 (2004), 365–392.
[5] B. Jüttler, C. Mäurer: Cubic Pythagorean-hodograph spline curves andapplications to sweep surface modeling. Computer-Aided Design, 31(1999), 73–83.
[6] S-H. Kwon: Solvability of G1 Hermite interpolation by spatial
Pythagorean-hodograph cubics and its selection scheme. Computer
Aided Geometric Design, 27 (2010), 138–149.
[7] B. O’Neill: Elementary Differential Geometry (2nd ed.), Academic
Press, 1997.
[8] F. Pelosi, R.T. Farouki, C. Manni, A. Sestini: Geometric Hermite interpolation by spatial Pythagorean-hodograph cubics. Advances in Computational Mathematics, 22 (2005), 325–352.
[9] M.G. Wagner, B. Ravani: Curves with rational Frenet-Serret motion.
Computer Aided Geometric Design, 15 (1997), 79–101.

138

Free vibration of symmetric angle-ply laminatedcylindrical shells of variable thickness including shear
deformation theory: spline method
K.K.Viswanathan*, Zainal Abdul Aziz, Saira Javed
Department of Mathematics, Faculty of Science, Universiti Teknologi Malaysia,
81310 Johor Bahru, Johor, Malaysia* e-mail: visu20@yahoo.com
Abstract
The vibrational study of structural elements is of persistent interest in theareas of missiles, aviation, shipping, surface transport and a large number ofindustries related to the cement and chemicals. Shells of layered construction
possess more desirable characteristics over those of homogeneous structure, interms of higher specific stiffness and better damping and shock absorbing
characteristics.Variation in thickness of the layers and inclusion of shear
deformation and rotary inertia provides further facilities for the designer inrealizing the desired vibrational characteristics of the structure. Transverse effects
become more pronounced as the shell become thicker relative to its in-plane
dimension and radius of curvature. Free vibrational study of symmetric angle-ply
laminated cylindrical shells of variable thickness including first order shear
deformation theory using spline function approximation is studied. The researchers[1-3] worked on cross-ply and angle-ply plates with constant thickness and conicalshells with variable thickness neglecting the shear deformation and rotary inertia.
The available works which consider also these complexities are very few and, tothe authors’ knowledge, do not use the versatile spline function method of solutionsand variable thickness [4-5]. In the present work, the results are expected to be
more accurate and interesting, and more suitable for immediate application in theareas already explained elsewhere.
The equations of motion for the cylindrical shells are derived using first order
shear deformation theory. The solutions of displacement functions are assumed in a
separable form to obtain a system of coupled differential equations in terms

139displacement and rotational functions and these functions are approximated by
Bickley-type splines of order three. The vibrations of three and five layered shells,
made up of several types of layered materials and subjected to two types ofboundary conditions are considered. A generalized eigenvalue problem is obtained
and solved numerically for an eigenfrequency parameter and an associatedeigenvector of spline coefficients. Parametric studies are made for the frequency
parameters with respect to the coefficients of thickness variations, length-to-radiusratio, length-to-thickness ratio, various types of materials and ply angles under
different boundary conditions.

References[1] K.K.Viswanathan, S.K. Lee: Int. J. Mechanical Science, 49 (2007),
352—363.
[2] K.K.Viswanathan, K.S. Kim: Int. J. Mechanical Sciences, 50 (2008),
1476—1485.
[3] K.K.Viswanathan, P.V. Navaneethakrishnan: Int. J. Solids & Structures,
42 (2005), 1129--1150.
[4] U. Topal, U. Uzman: Materials & Design, Vol.30 (2009), 532--536.
[5] T.Timarci, K.P. Soldatos: J. Engineering Mathematics, Vol.37 (2000),
211--230.

140

Point vortex equilibria enhancing forces overtwo parallel plates
Takashi Sakajo1),∗
1) Department of mathematics, Hokkaido University,
CREST, Japan Science and Technology Agency∗ e-mail: sakajo@math.sci.hokudai.ac.jp
Abstract
In ’70s, Kasper designed a high aspect ratio wing with a ﬂap, which
is known as “Kasper airfoil”. He suggested that the lift of the airfoil wasenhanced by placing vortex structures around its wing. The concept has been
examined by some wind-tunnel experiments[3, 7]. Rossow[4] demonstrated
that it was possible to keep a spanwise vortex line coupled with an axial ﬂow
along its core over the wing with a vertical ﬂap. Saﬀman and Sheﬃeld[5]considered a potential ﬂow model describing a free point vortex over a ﬂatplate in a uniform ﬂow. They found that some of stationary point vortices
can stay stably above the plate and they enhance the lift without any ﬂap.
In the present talk, we consider a ﬂow over two parallel ﬂat plates in thepresence of a uniform ﬂow whose schematic picture is shown in Figure 1. Thisis an extension of the analysis due to Saﬀman and Sheﬃeld into a doubly
connected ﬂow domain to ﬁnd vortex equilibria enhancing the lift on the
double wing. It is also motivated by a study of an design for vertical axis
wind turbines, which are wind-powered electronic generators with verticalmultiple blades[8]. With keeping the same concept as the Kasper airfoil
in mind, we expect that stationary vertical vortex lines could enhance therotational force when they were attached stably in the neighborhood of theblades. As a simple two-dimensional model, we consider two stationary point
vortices around the two parallel plates and observe how the stationary point
vortices enhance the forces acting on the two plates.
The analysis is based on a calculus for potential ﬂows in multiply connected domains introduced by Crowdy[1]. The equation of motion for two
point vortices is derived from a mathematical formulation of the N -vortex
problem in multiply connected domains given in [6]. We ﬁnd that there exist
loci of stationary two point vortices that enhance the upward lift and theclockwise rotation of the two parallel plates.

141

S
ʣ

w２

ʣ

u１

0q

Figure 1: We consider ﬁxed conﬁgurations of two point vortices behind twoparallel plates in a uniform ﬂow and observe how the ﬁxed equilibria enhance
the force on the plates.

References[1] D. G. Crowdy, A new calculus for two dimensional vortex dynamics,
Theor. Comput. Fluid Dyn., 24 9–24 (2010) doi:10.1007/s00162-0090098-5[2] D. Hummel, On the vortex formation over a slender wing at large angles
of incidence, AGARD CP-247, paper no. 15.
[3] E. W. Kruppa, A wind tunnel investigation of the Kasper vortex concept,
AIAA paper, 77-310, Jan. (1977)[4] V.G. Rossow, Lift enhancement by an externally trapped vortex, AIAA
paper 77-672, June (1977)[5] P. G. Saﬀman and J. S. Sheﬃeld, Flow over a wing with an attachedfree vortex, Stud. Appl. Math., 57 107–117 (1977)[6] T. Sakajo, The N -vortex problem in multiply connected domains, Proc.
Roy. Soc. A, 465 2589–2611 (2009)[7] D. Walton, A brief wind tunnel test of the Kasper aerofoil, Soaring, 38
No.11 26–27 (1974)[8] Seethewebpageat
http://en.wikipedia.org/wiki/Verticalaxis wind turbine for a description of vertical axis wind turbines.

142

Bifurcation analysis for the Lugiato-Lefever
equation in a square
Tomoyuki Miyaji1),∗ , Isamu Ohnishi2) , Yoshio Tsutsumi3)
1) Research Institute for Mathematical Science, Kyoto University2) Department of Mathematical and Life Sciences, Hiroshima University3) Department of Mathematics, Kyoto University∗ e-mail: tmiyaji@kurims.kyoto-u.ac.jp
Abstract
We show some results on bifurcation analysis and numerical simulation
for the Lugiato-Lefever equation (LLE) on a square domain. LLE is a cubic
nonlinear Schrödinger equation (NLS) with damping, detuning and driving
force. It is a model for describing the evolution of transversal patterns in an
optical cavity with a Kerr medium [1]. It is given by∂t E = − (1 + iθ) E + ib2 ∆E + i|E|2 E + Ein ,x ∈ Ω, t > 0,(1)

where x = (x1 , x2 )T , Ω = (−1/2, 1/2)2 ⊂ R is a unit square, and ∆ =∂x21 + ∂x22 is the Laplacian. We impose the periodic boundary conditions. The
parameters b2 , θ ∈ R are diﬀraction and detuning parameters, respectively.
E = E(x, t) ∈ C denotes a slowly varying envelope of electric ﬁeld. Ein ≥ 0denotes the intensity of spatially homogeneous driving ﬁeld, and it is a main
control parameter.
Numerical simulations suggest that LLE in 1- or 2-dimensional space has
solitary wave solutions in a certain range of parameters[2]. In contrast to
NLS, LLE does not satisfy any conservation law of NLS. Typically, such asolution of LLE appears as an equilibrium point. Moreover, for 2D LLE, alocalized spot can undergo a Hopf bifurcation, and it results in a spatiallylocalized and temporally oscillating solution called oscillon [3]. Such solutions
are understood as dissipative structures resulting from the balance between
gain and loss of energy. However, no explicit analytical solutions are known
for LLE. The authors have studied steady-state bifurcations for 1D LLE [4].
They have proved that a small localized roll can bifurcate from homogeneousstate in R, and that two mixed-mode solutions can bifurcate as a secondary
bifurcation in T1 , which are considered to be “germs” of localized structures.
We are interested in oscillons for 2D LLE. As a ﬁrst step, we study steadystate bifurcation of spatially homogeneous steady state in a mathematicallyrigorous sense, because it is responsible for the occurrence of stationary localized patterns. For this purpose, we apply the center manifold reduction
and group theoretic bifurcation theory [5].
143We study (1) with periodic boundary conditions. Taking into account of
the symmetry of the system, we can obtain a catalog of model-independent
bifurcation behaviors. Then we apply the center manifold reduction near
the modulational instability of homogeneous state. The symmetry of LLE
restricts the function form of the vector ﬁeld on the center manifold. It helps
to study model-speciﬁc behaviors. The symmetry of LLE is described by
Γ = D4 n T2 , where n means semidirect product of groups. D4 is generated
by the reﬂection across x2 = x1 and π/2-rotation about the origin. T2 is agroup of translations modulo 1.
Suppose that k = (l, n) ∈ Z2 is a critical wave vector of modulational
instability. We may assume l ≥ n ≥ 0 without loss of generality. Since weconsider single-mode bifurcations, we assume that the solution (l, n) ∈ Z2
to l2 + n2 = k for a given k ∈ N is unique in the above sense, There arethree cases to be distinguished: i) l = n > 0; ii) l > n = 0; iii) l > n > 0.
In the ﬁrst two cases, the center manifold is 4-dimensional, while that is8-dimensional in the third case. We classify the possible bifurcations for(1) in a small neighborhood of the bifurcation point in a mathematicallyrigorous sense. Especially, it turns out that all bifurcating solutions are
unstable near the bifurcation point. However, numerical simulations suggestthat there exist stable patterns. In order to capture them, we have to study
bifurcations with higher codimension.

References[1] L. A. Lugiato and R. Lefever, Spatial Dissipative Structures in Passive
Optical Systems, Phys. Rev. Lett., 58 (1987), 2209–2211.
[2] A. J. Scroggie, W. J. Firth, G. S. McDonald, M. Tlidi, R. Lefever and L.
A. Lugiato, Pattern formation in a passive kerr cavity, Chaos, Solitons& Fractals, 4 (1994), 1323–1354.
[3] P. Colet, D. Gomila, A. Jacobo and M. A. Matı́a, Excitability mediated by dissipative solitons in nonlinear optical cavities, Lect. Notes inPhysics, 751 (2008), 113–135.
[4] T. Miyaji, I. Ohnishi and Y. Tsutsumi, Bifurcation analysis to the
Lugiato-Lefever equation in one space dimension, Physica D, 239 (2010),
2066–2083.
[5] M. Golubitsky, I. Stewart and D. G. Schaeﬀer, Singularities and Groups
in Bifurcation Theory II, Applied Mathematical Sciences, 69, SpringerVerlag, New York (1988).
144

A theorem on S-shaped bifurcation curve fora positone problem with convex-concave
nonlinearity and its applications to theperturbed Gelfand problem
Kuo-Chih Hung*, Shin-Hwa Wang
Department of Mathematics, National Tsing Hua University,
Hsinchu, Taiwan∗ e-mail: kchung@mx.nthu.edu.tw
Abstract
We study the bifurcation curve and exact multiplicity of positive solutionsof the positone problem½ 00
 () + () = 0, − 1    1,
(−1) = (1) = 0,
where   0 is a bifurcation parameter,  ∈  2 [0 ∞) satisfies  (0)  0 and()  0 for   0, and  is convex-concave on (0 ∞) and is asymptotic
sublinear. Under a mild condition, we prove that the bifurcation curve is
S-shaped on the ( kk∞ )-plane. We give an application to the perturbed
Gelfand problem¡  ¢½ 00= 0, − 1    1,
 () +  exp +(−1) = (1) = 0,
where   0 is the Frank—Kamenetskii parameter and   0 is the activationenergy parameter. We show that, if  ≥ ∗ ≈ 4166 for some constant ∗ , thebifurcation curve is S-shaped on the ( kk∞ )-plane. Our results improvethose in Korman and Li [1] and Wang [2].

References[1] P. Korman, Y. Li, On the exactness of an S-shaped bifurcation curve,
Proc. Amer. Math. Soc. 127 (1999) 1011—1020.
[2] S.-H. Wang, On -shaped bifurcation curves, Nonlinear Anal. 22 (1994)
1475—1485.

145

On the existence of a double S-shapedbifurcation curve with six solutions
for a combustion problem
Kuo-Chih Hung, Shin-Hwa Wang , Chien-Shang Yu
Department of Mathematics, National Tsing Hua University,
Hsinchu, Taiwan
e-mail: shwang@math.nthu.edu.tw
Abstract
We study the combustion problem with nonlinear boundary conditions
given by8u00>
< u (x) = exp( +u ), 0 < x < 1,u(0) = 0,
(1)
>
u(1)
: u(1) u0 (1) + [1
]u(1)
=
0,
u(1)+1u(1)+1where> 0 is the Frank–Kamenetskii parameter or ignition parameter,a > 0 is the activation energy parameter, u(x) is the dimensionless temperau) is the temperature dependence obeying
ture, and the reaction term exp( +uthe simple Arrhenius reaction-rate law in irreversible chemical reaction kinetics. We mainly prove that, for > 1 6:459 for some constant 1 , thebifurcation curve of the problem is double S-shaped on the ( ; kuk1 )-plane
and the problem has at least six positive solutions for a certain range of .
We give rigorous proofs of some computational results of Goddard II, Shivaji
and Lee [1].

Fig. 1. Double S-shaped bifurcation curve of (1) for

146

>

1

6:459.

Theorem 1. (See Fig. 1.) Consider problem (1). Suppose
>
6:459
for
some
constant,
thenthebifurcation
curveofproblem(1)1
1
is double S-shaped on the ( ; kuk1 )-plane. Moreover, there exist positive
numbers 3 < 1 < 0 < 4 < 2 such that(i) If0 <(ii) If

1 <(iii) If

4 <(iv) If(v) If<

4 , problem (1) has at least 6 positive solutions,
0 or=

4 , problem (1) has at least 5 positive solutions,

<

2 or

=

1 , problem (1) has at least 4 positive solutions,

3 <

<

1 or

=

2 , problem (1) has at least 3 positive solutions,

>

2 or

=

(v) If 0 <<3 , problem (1) has at least 2 positive solutions,

3 , problem (1) has a unique positive solutions.

References[1] J. Goddard II, R. Shivaji, E. K. Lee: A double S-shaped bifurcation
curve for a reaction-di¤usion model with nonlinear boundary conditions,
Boundary Value Problems, 2010, Article ID 357542, 23 pages.

147

Appearance and Disappearance ofNon-symmetric Progressive Waves of Deep Water
Mayumi SHŌJI
1),2)

1),*

2)

, Chika SHIMIZU

Department of Mathematical and Physical Sciences
Japan Women’s University* e-mail: shoji@fc.jwu.ac.jp
Abstract

We consider a problem of progressive water waves with a constant speed. It is
supposed to be two-dimensional and irrotational flow of incompressible inviscid
fluid. Under the effect of gravity and surface tension, it becomes the bifurcation
problem which has two bifurcation parameters concerned with these two forces
respectively. It is a classical problem and there have been a lot of researches
mathematically and numerically. However most of them are restricted to symmetricwaves and we know a few about non-symmetric waves. Numerically there are only
the works of Zufiria for gravity waves on water of finite depth [2], infinite depth[3] and capillary-gravity waves of finite depth [4]. All of them are the waves with
six peaks in one wave length. We made a thorough investigation in case ofcapillary-gravity waves of infinite depth. We confirmed the existence of similar
non-symmetric waves to [4] and saw the mechanism of symmetry breaking. And
furthermore we were able to obtain some other new types of non-symmetric waves.
       
We studied the problem using the formulation with Levi-Civita’s equation. And
numerically we applied the spectral-collocation method and Keller’s method (see[1] for details). The formulation and the numerical methods we used here are quite
(A)

(B)

Figure 1. p=1.5. (A) The dotted line is the branch of non-symmetric waves. The
branches of symmetric waves are just on the q-axis. (B) The wave profile of asolution at (b) on the branch of non-symmetric waves.
148different from [2-4]. Figure 1 shows an example of our numerical results of
non-symmetric waves with six peaks. The left figure is the bifurcation diagram,where the vertical axis q is the bifurcation parameter related to surface tension. The
horizontal axis b2 represents the asymmetric degree of non-symmetric waves,namely b2=0 if wave is symmetric. Another bifurcation parameter p related togravity is fixed here. The branch of non-symmetric waves constitutes a closed loop
(dotted line in (A)) and the branches of symmetric waves are just on the q-axis. The
symmetry-breaking occurs in the following mechanism: a trivial solution  asimple bifurcation of mode 6 symmetric waves  a secondary bifurcation ofmode 2  a tertiary bifurcation of mode 1  symmetry-breaking. We can see
the looped branch of non-symmetric waves becomes smaller as p decreases anddisappears when p  1.3.
     
Other than the above type of waves with six peaks, we were able to obtain some
new types of non-symmetric waves. One is a type of non-symmetric wave withthree peaks in one wave length and another one is a type with seven peaks. Maybe,
if any other non-symmetric solution exist, we predict that might be the bifurcation
solution of a higher mode than six at first. Hence it is unexpected and interestingfor us to obtain non-symmetric waves of mode 3.

Figure 2. Non-symmetric waves of mode 3 and mode 7.

References[1] H. Okamoto and M Shōji, The Mathematical Theory of Permanent
Progressive Water-Waves, World Scientific (2001).
[2] J. A. Zufiria, Weakly nonlinear non-symmetric gravity waves on water offinite depth, J. Fluid Mech., 180 (1987), pp. 371--385.
[3] J. A. Zufiria, Non-symmetric gravity waves on water of infinite depth, J.
Fluid Mech., 181 (1987), pp. 17--39.
[4] J. A. Zufiria, Symmetry breaking in periodic and solitary gravity-capillary
waves on water of finite depth, J. Fluid Mech. 191 (1987), pp. 341--372.

149

Bifurcation structure and spontaneous
pattern formation in binary fluid convection
Takeshi WATANABE1),∗ , Makoto IIMA1),2) , Yasumasa NISHIURA1)
1) Research Institute for Electronic Science, Hokkaido University
2) JST, CREST∗ e-mail: watanabe@nsc.es.hokudai.ac.jp
Abstract
We study dynamical behavior of localized pattern in a binary ﬂuid convection. In a binary ﬂuid mixture, thermal convection exhibits localized patterns
such as steady pulse (SP), a stationary spatially localized convection cell; andperiodic traveling pulse (PTP), a moving spatially localized convection cell
whose shape periodically changes with time, as well as ordinary RayleighBénard type convection. So far, PTP solution has not been obtained evennumerically because the solution has two unknown variables: group velocity
and temporal period in the comoving frame with the group velocity. It is
important to obtain such solutions including their stability in both views ofdetermining mathematical structure of this system and controlling numericalexperiments. We succeeded in obtaining bifurcation structure of PTP, SP,
and SOC solutions[2].
The method to obtain PTP solutions is established by setting an appropriate dimensional Poincaré section. For N -dimensional dynamical system,
PTP solution can be determined uniquely by setting (N − 2)-dimensional
Poincaré section because PTP solution have time and spatial translational
invariance. In other words, a point on a torus spanned by PTP solution vector is identiﬁed uniquely by two values. Therefore we search for a solution to
the N -dimensional dynamical system in (N − 2)-dimensional space and thus
the number of unknowns does not equal the number of equations, however,this lack is compensated by two additional unknowns, i.e., group velocity andperiod, and then a closed set of equations is obtained. The Newton-Raphson
method is used to ﬁnd the ﬁxed point on the Poincaré section together with
the GMRes method[2].
The obtained global bifurcation diagram for Γ = 32 revealed that bifurcation branch of SP forms a snake-like structure which is frequently observed inmany dissipative systems including a plane Couette ﬂow[4] and a reactiondiﬀusion system[5]. Thus, various SP solutions with diﬀerent numbers ofconvection cells coexist in a parameter range. Further, PTP solutions also
coexists in the same parameter range. This fact is consistent with the result150

of collision analysis of two counter-propagating PTPs reported by Toyabe[3],in which various states are obtained as the result of the collision depending on the phase of each PTPs. However, he obtained PTP by calculating
asymptotic state of time evolution, that is, the accuracy was not enough toinvestigate the detailed input-output relationships of the collision. In thepresent study, we investigate detailed phase dependency of the PTP-PTP
collision by controlling obtained solutions. Surprisingly, almost all cases result in a single PTP in the asymptotic state and SP states are observed fora few cases where the initial value is close to mirror-symmetric. This result
supports the robustness of PTP.
We also study how PTP is generated. Since there exist stable SP and PTP
branches even though the conductive state is unstable, it is suggested thatexistence and stability of these solutions depend on the size of the system.
We performed time evolution in a very large domain (Γ = 500) to clarify thesize eﬀect and found that PTP and SP lose their stability even if they arestable when Γ = 32. The convective instability contribute to the perturbation
growth, which is consistent with the result of Batiste et al[1]. After complex
interactions, we observed several PTPs propagating in the same direction
with non-uniform intervals.
A part of this study is supported by Grant-in-Aid for Scientiﬁc Research(KAKENHI) No. 21340019 and Core Research for Evolutional Science andTechnology (CREST) No. PJ74100011.

References[1] O. Batiste, E. Knobloch, A. Alonso, and I. Mercader, Spatially localizedbinary-ﬂuid convection, J. Fluid Mech. 560, 149–158 (2006).
[2] T. Watanabe, K. Toyabe, M. Iima, and Y. Nishiura, Time-periodic traveling solutions of localized convection cells in binary ﬂuid mixture, Theoretical and Applied Mechanics Japan, 59, 211–219 (2010).
[3] Kazutaka Toyabe, Collision dynamics of localized convection cells inbinary fluid mixture: network structure of collision orbit, Master Thesis,
Graduate School of Science, Hokkaido University (2009) [in Japanese].
[4] T. M. Schneider, J. F. Gibson, and John Burke, Snakes and ladders:localized solutions of plane Couette ﬂow, Phys. Rev. Lett. 104, 104501(2010).
[5] Y. Nishiura, D. Ueyama, Spatio-temporal chaos for the Gray-Scottmodel, Physica D, 150, 137–162 (2001).
151

Single-machine group scheduling problems with
time-dependent learning effect and position-based setup
time learning effectWen-Hung Kuo
Department of Information Management, National Formosa University
e-mail: whkuo@nfu.edu.tw
Abstract
Learning and its effect on productivity is well recognized in industry. In general,
the worker learns from experience in performing repetitive tasks and the learning
effect manifests itself in labor productivity improvement. Therefore, schedulingresearch has increasingly taken the concept of learning into consideration. There
are two main classes of learning effect models studied in the scheduling literature.
One is position-based learning. For example, Biskup [1] introduced a schedulingmodel with a job-independent learning effect in which the actual processing time of
job J i is pir = pi r a if it is scheduled in position r of a sequence, where pi anda ( a < 0 ) are the normal processing time and the learning rate of job J i ,respectively. He showed that the single-machine scheduling problems with thelearning effect to minimize the deviation from a common due date and the total
completion time are polynomially solvable. Mosheiov and Sidney [2] considered amore general model of learning effect in which the learning rates of the jobs aredifferent, i.e., the learning effects are job-dependent so that pir = pi r ai . They
showed that some scheduling problems with the job-dependent learning effect arepolynomially solvable.
The other is sum-of-processing-time-based learning. For example, Kuo andYang [3] introduced a scheduling model with a time-dependent learning effect ina r −1which the actual processing time of job J i is pir = pi 1 + ∑ p[ k ]  if it is k =1scheduled in position r of a sequence, where pi and p[ k ] are the normal

processing times of job J i and the job scheduled in position k, respectively, and ais the learning rate. They showed that the problem to minimize total completiontime is polynomially solvable. Koulamas and Kyparisis [4] proposed anothersum-of-processing-time-based learning model in which the actual processing timeα ∑ r −1 p[ k ] of job J i is pir = pi 1 − kn =1
 ∑ p[ k ] 
k =1



if it is scheduled in position r of a sequence,152where pi and p[ k ] are the normal processing times of job J i and the job
scheduled in position k, respectively, and α ≥ 1 . They also provided polynomial
solutions for some scheduling problems. For a comprehensive survey on thisstream of scheduling research, the reader may refer to Biskup [5].
Biskup [5] stated that it may be more appropriate to assume that position-basedlearning
takes
place
duringmachinesetupsonly,
while
sum-of-processing-time-based learning occurs in considering the experience thatworkers have gained from producing jobs.
To the best of our knowledge, there is still no study focusing on the viewpointin the survey paper [5]. Therefore, in this paper, we will consider
sum-of-processing-time-based learning on job processing time and position-basedlearning on setup time in some single-machine group scheduling problems. The
objective is to minimize the makespan. We provide two polynomial timealgorithms to solve the makespan minimization problems.

References[1] D. Biskup. Single-machine scheduling with learning considerations,
European Journal of Operational Research, 115 (1999), 173-178.
[2] G. Mosheiov, J.B. Sidney. Scheduling with general job-dependent learning
curves, European Journal of Operational Research, 147 (2003), 665-670.
[3] W.H. Kuo, D.L. Yang. Minimizing the total completion time in asingle-machine scheduling problem with a time-dependent learning effect,
European Journal of Operational Research, 174 (2006), 1184-1190.
[4] C. Koulamas, G.J. Kyparisis. Single-machine and two-machine flowshop
scheduling with general learning functions, European Journal of Operational
Research, 178 (2007), 402-407.
[5] D. Biskup. A state-of-the-art review on scheduling with learning effects.
European Journal of Operational Research, 188 (2008), 315-329.


153

Optimal Retirement with Unemployment
Risks and Market Completion
Bong-Gyu Jang∗ , Seyoung Park, Yuna Rhee
Department of Industrial and Management Engineering, POSTECH∗ e-mail: bonggyujang@postech.ac.kr
Abstract
This paper investigates individual’s optimal retirement in the presenceof involuntary unemployment risks and positive wealth constraint. Most
existing literature on optimal retirement problem, e.g., the paper of Farhi
and Panageas [1], tries to find individual’s optimal consumption, investment,
and voluntary retirement time simultaneously. We add the possibility ofher involuntary unemployment to the conventional set-up, and assume thatthe involuntary unemployment readily breaks out when a crucial exogenous
shock such as the bankruptcy of the company for which she is working comes
into being. To reach the goal, we introduce a complete market with frictions.
In this regard we consider a complete market with a personalized unemployment insurance, and accordingly, the individual can hedge the involuntary
unemployment risks by purchasing the insurance policy. However, we assume
that the price of the insurance policy is costly compared with its actuariallyfair price to reflect the current situation that the unemployment insurance
markets are quite thin.
We utilize the martingale approach to take into account an optimal stopping for voluntary retirement under the conditions above. The technicalcontribution of our paper is in this point of view; to find an analytic solution of the optimal retirement problem in the presence of the uncertainties ofretirement time by using the martingale approach. In our model, since thestopping time of the optimal retirement should be determined endogenously
under the condition that the individual might lose her job accidentally, one
can say that our paper is an extension of the work of Farhi and Panageas[1]. We indeed do such job even with individual’s positive wealth constraint,
which makes our model go near the real world and precludes arbitrage opportunities. Even though a plethora of studies explored individual’s optimal
retirement, as far as we know, this is the first study examined not only withthe voluntary retirement but also with the involuntary unemployment. More
importantly, it allows one to capture more various economic features thanthose of Farhi and Panageas [1].
Basically, in our model, a wealth process equipped with an unemployment
insurance is considered. When πt , ct , I, k, τ , τU , and W t present the dollar154amount of investment in the stock, consumption rate, income rate, fixed rateof the difference between the market price and actuarially fair price of theunemployment insurance rate, voluntary retirement time, involuntary unemployment time which is governed by exponential distribution with intensity
λU , and potential voluntary retirement wealth level, respectively, the wealthprocess Wt follows byhidWt = rWt + πt (µ − r) − ct + I1{t<τ ∧τU } − k − λU {W t − Wt } dt + πt σdBt ,where r and σ denotes the mean rate of return and volatility of a stock,
respectively. Note that, to reflect the current unemployment insurance markets, we introduce the additional rate k representing market frictions in trading the unemployment insurance policy which can leads to the effect of reduction in income I. The insurance policy helps to achieve a complete market inthat involuntary unemployment risks can be hedged through the instrument.
Besides, an individual gets unemployment coverage purchased for a premiumof λU (W t − Wt ) per unit time and a receipt of W t − Wt at unemployment.
We find some interesting properties are related to λU and k in the following view points: (1) the optimal retirement wealth level, (2) optimal
consumption, (3) optimal risky portfolio, (4) the expected time to retirement, and (5) the probability that the individual enters voluntary retirement
before involuntary unemployment. The high intensity of involuntary unemployment and high cost of the insurance policy tend to lead the individual
tries to enter voluntary retirement earlier, which causes the lower retirement
wealth level and less amount of consumption but the larger proportion ofstockholdings. This is because she gets lower income due to the cost andwants to hedge against the involuntary unemployment risks. Moreover, theexpected time to retirement decreases as the wealth level increases since thewealth level increasingly close to threshold wealth level. Thus, the higher
intensity of involuntary unemployment and cost of the insurance policy, theshorter expected time to retirement. Finally, growing wealth level makes theprobability that the individual enters voluntary retirement before involuntary unemployment high, and the probability rapidly increases as the highintensity and cost.

References[1] E. Farhi, S. Panageas: Saving and Investing for Early Retirement: A
Theretical Analysis, Journal of Financial Economics, 83 (2007), 87–121.

155

Rescaled Gradient-based Methods with Fixed
Stepsize for Large-scale Optimization
Wah June Leong1,2),∗ , Mahboubeh Farid2) , Malik Abu Hassan1)
1) Department of Mathematics, Faculty of Science, University Putra
Malaysia
2) Institute for Mathematical Research, University Putra Malaysia∗ e-mail: leong@math.upm.edu.my
Abstract
A new family of gradient-type method with improved Hessian approximation for large-scale unconstrained optimization problems is introduced. The
new method resembles the Barzilai-Borwein (BB) [1] method, except thatthe Hessian matrix is approximated by a diagonal matrix rather than themultiple of the identity matrix in the BB method. Then the diagonal Hessian approximation is derived based upon the commonly used least change
updating strategy for quasi-Newton methods with the added restriction thatfull matrices are replaced by diagonal matrices [2]. Using some appropriate
matrix norms, some members of this family are introduced. Numerical experiments with some particular members of the family are performed andthe results show that the proposed method yields desirable improvement.

References[1] J. Barzilai and J.M. Borwein: Two point step size gradient methods,
IMA J. Numer. Anal. 8 (1988) 141-148.
[2] M.A. Hassan, W.J. Leong and M. Farid, A new gradient method
via quasi-Cauchy relation which guarantees descent, J. Comput. Appl.
Math. 230 (2009) 300-305.

156

Behavior of polygonal curves by crystalline
curvature ﬂow with a driving force
Tetsuya Ishiwata1),∗
1) Department of Mathematical Science, Faculty of System Science andTechnology, Shibaura Institute of Technology∗ e-mail: tisiwata@shibaura-it.ac.jp
Abstract
In this presentation we consider the motion of planar polygonal curves
governed by generalized crystalline curvature ﬂow with a driving force:β(Nj )Vj = U − g(Hj ),(1)

where Vj , Nj and Hj denote an outward velocity, an outward normal vector
and a crystalline curvature of the j-th facet of solution curve, respectively.
Here, “facet” means a lines segment of solution curve Γ(t) and Fj denotes
the j-th facet of Γ(t). The positive function β and U describe an anisotropyof the mobility and a driving force, respectively. In this paper we consider
the case that U > 0 and the function g = g(λ) is odd, locally Lipschitz
continuous, monotone increasing in λ and has linear or superlinear growthfor λ ∼ ±∞.
The equation (1) with g(λ) = λ is a simple model of an interface motionof crystals (See [1] and [6]). In the case when the solution region is bounded,
many authors discuss the behavior of the solutions for several types of crystalline motions (See [2, 7, 3] and their references). On the other hands, in thecase when the solution region is unbounded, there are a few results. Marutani et al. [5] shown the existence of traveling “V-shaped” solutions for thecase g(λ) = λ as the limit of the smooth solution for weighted curvature ﬂow.
In [4], we consider the simple case β ≡ 1 and g(λ) = λ and show that thesolution with non-V-shaped initial region eventually become V-shaped undersome conditions. This kind of monotonicity phenomena are already known
for the case when the solutions are bounded. The most famous result isshown by M. Grayson in 1987 for smooth solution curves to curve-shortening
ﬂow. He showed that the solution curves from non-convex initial curve becomes convex in ﬁnite time. For crystalline case, it is shown in [3] for thecase U = 0 that crystalline curvature of all facets of the solution polygons
from non-convex initial polygon becomes non-negative in ﬁnite time.
In this presentation we show the extended results in [3, 4] for generalized
motion (1).
157

References[1] S. Angenent and M.E. Gurtin: Multiphase thermomechanics with interfacial structure, 2. Evolution of an isothermal interface, Arch. Rational
Mech. Anal. 108 (1989), 323–391.
[2] Y. Giga and M. E. Gurtin: A comparison theorem for crystalline evolution in the plane, Quart. J. Appl. Math. LIV (1996), 727–737.
[3] T. Ishiwata: Motion of non-convex polygons by crystalline curvature andalmost convexity phenomena, Japan Journal of Industrial and Applied
Mathematics, 25 (2008), no. 2, 233–253.
[4] T. Ishiwata: On the motion of polygonal curves with asymptotic lines
by crystalline curvature ﬂow with bulk eﬀect, to appear in Discrete andContinuous Dynamical Systems, Series S.
[5] Y. Marutani, H. Ninomiya and R. Weidenfeld: Traveling curved frontsof anisotropic curvature ﬂows, Japan Journal of Industrial and Applied
Mathematics, 23 (2006), no. 1, 83–104.
[6] J. E. Taylor: Constructions and conjectures in crystalline nondiﬀerential
geometry, Proceedings of the Conference on Diﬀerential Geometry, Riode Janeiro, Pitman Monographs Surveys Pure Appl. Math. 52 (1991),
321–336, Pitman London.
[7] S. Yazaki: Point-extinction and geometric expansion of solutions to a
crystalline motion, Hokkaido Math. J. 30 (2001), 327–357.

158

Some remarks on the fictitious domain
method with penalty for elliptic problems
Guanyu ZHOU1),∗ , Norikazu SAITO1)
1) Graduate School of Mathematical Sciences, The University of Tokyo∗ e-mail: zhoug@ms.u-tokyo.ac.jp
Abstract
The principle of the fictitious domain method is to solve the problem in alarger domain (the fictitious domain) containing the domain of interest witha very simple shape. Then, the fictitious domain is discretized by a uniform
mesh, independent of the original boundary. The advantage of this approachis that we can avoid the time-consuming construction of a boundary-fitted
mesh. One of these approaches is the penalty fictitious domain method
which is based on a reformulation of the original problem in the fictitious
domain by using penalty parameter. Obviously, this approach is of use totreat time-dependent moving-boundary problems. However, it seems thatpenalty fictitious domain method is not studied very well at present. The
main purpose of this work is to study the penalty fictitious domain method
applied to these time-dependent moving-boundary problems. As a primarystep towards this final end, herein we examine some new methods of error
analysis for elliptic problems that can be applied to time-dependent problems
with no much difficulties.
We consider the original elliptic problem (Q),
Find u ∈ H01 (Ω) such that(1)(∇u, ∇v)Ω = (f, v)Ω , ∀v ∈ H01 (Ω),
where Ω ⊂ R2 denotes a smooth bounded domain, (·, ·)Ω is the inner productof L2 (Ω) and f ∈ L2 (Ω). We can find a rectangular domain D ⊃ Ω, Ω1 =D\Ω, and turn to solve the approach problem (Qǫ ) with penalty coefficient
0 < ǫ ≪ 1,
Find uǫ ∈ H01 (D) such that(2)(∇uǫ , ∇v)Ω + 1ǫ (∇uǫ , ∇v)Ω1 = (f˜, v)D , ∀v ∈ H01 (D),where f˜ is the zero extension of f into D.
The error kuǫ − uk1,Ω is analyzed by many√authors, where k · k1,Ω denotes
the H 1 (Ω) norm. In [1], it is bounded by C ǫ (C is some constant, so as
in the following), and in [2, 3] the sharper estimate Cǫ is achieved. In this159

paper, we give a new way to derive the sharper estimates, as well as someother estimates, such as kuǫ k1,Ω1 ≤ Cǫ. We find our method of analysis iseasy to be applied to parabolic problems. Moreover, we present some otheranalysis and estimates of uǫ , which is useful for studying the H 1 and L2
error between the solutions of approach problem (Qǫ ) and its discretization
problem, which is denoted as (Qǫ,h ).
A Cartesian mesh can be introduced to the rectangular domain D to get
a uniformed triangulation Th , h is the maximum diameter of the triangles ofTh . Vh (D) is a subspaces of piecewise linear continuous functions subordinate
to Th . Then (Qǫ,h ) reads as:Find uǫ,h ∈ Vh (D) such that(3)(∇uǫ,h , ∇vh )Ω + 1ǫ (∇uǫ,h , ∇vh )Ω1 = (f˜, vh )D , ∀vh ∈ Vh (D).
In the literature, there are numerous works devoted to the study of the H 1
error between uǫ,h and uǫ in√Ω. For example,in [2], it is proved that the H 1√ √error is bounded by C(ǫ + ǫ h + h ǫ + h). In our work, we prove√ a similar
result. Moreover, we show the L2 error is bounded by C(ǫ + h + ǫh).
In above, we have restricted our attention to Dirichlet boundary problem.
For Neumann and mixed boundary problems we also study the approximation of approach problems( they are a little different from the Dirichlet case),
as well as the discretization problem for Neumann case. Although the main
results we obtain are much similar to those of [2], we omit to present thedetail here. However, as we mentioned in the beginning, our work is focus on solving parabolic problem with time-dependent domain, and all our
analyzing methods are applicable well.

References[1] A. L. Mignot: Méthodes d’approximation des solutions de certains
problèmes aux limites linéaires. I, Rend. Sem. Mat. Univ. Padova, 40(1968) 1-138.
[2] S. Zhang: Analysis of finite element domain embedding methods forcurved domains using uniform grids, SIAM J. Numer. Anal. 46 (2008),
no. 6, 2843-2866.
[3] B. Maury: Numerical analysis of a finite element/volume penalty
method, SIAM J. Numer. Anal. 47 (2009), no. 2, 1126-1148.

160

Analysis of the finite volume approximation
for a degenerate parabolic equation
Norikazu SAITO1),∗
1) Graduate School of Mathematical Sciences, The University of Tokyo∗ e-mail: norikazu@ms.u-tokyo.ac.jp
Abstract
We consider the initial-boundary value problem for a degenerate parabolic
equation, ut − ∆f (u) = 0 in Ω × (0, T ),

u = 0 on ∂Ω × (0, T ),u|t=0 = u0 (x) on Ω,(1)

where Ω ⊂ R2 denotes a polygonal domain, T an arbitrary positive constant,and f a non-decreasing continuous function deﬁned on R satisfying f (0) =0. We also consider the boundary value problem for a degenerate elliptic
equation,{u − λ∆f (u) = g in Ω,(2)
u = 0 on ∂Ω,where λ is a positive constant and g a function deﬁned on R.
As is well-known, Problem (1) describes, for instance, the ﬂow of homogeneous ﬂuid in porous media if
f (u) = u |u|γ−1(3)with γ > 1, the fast diﬀusion if (3) with 0 < γ < 1, and the two phase Stefan
problem in enthalpy formulation ifα(u + 1) (u ≤ −1)
f (u) = 0(−1 < u < 1)β(u − 1) (u ≥ 1)with α, β > 0.
The purpose of this paper is to report some operator theoretical propertiesof the ﬁnite volume method (FVM) in the sense of [3] applied to (1) and(2). As is well-known, L1 theory of Brezis and Strauss [1] is of great use todeal with this problem. In this paper, we shall see that FVM is a suitablediscretization method for (1) and (2) in the sense that the discrete version
161

of [1] can be applied. Consequently, we immediately deduce the generationof the nonlinear semigroup, namely, the unique existence of a solution to asemidiscrete (in space) FVM for (1). Then, we readily obtain stability resultsin L1 and L∞ , and order-preserving property for numerical solutions via thenonlinear semigroup theory. This is totally new approach to study FVM fordegenerate elliptic and parabolic problems.
FVM is a discretization method based on local conservation properties ofequations so that it is well suited for PDEs of conservation laws. Although
the range of application seems to be smaller than that of the ﬁnite element
method (FEM), FVM has its own advantages. For example, FVM naturally
satisﬁes the discrete maximum principle, if a linear diﬀusion is considered.
We recall that the discrete maximum principle in FEM holds only when someshape conditions on the triangulation are satisﬁed, and such a restriction often causes some diﬃculties. In this paper, we shall reveal another advantage
of FVM through the degenerate diﬀusion problems and the nonlinear semigroup theory.
Though we shall restrict our consideration to the two dimensional polygon
in this paper, it is not diﬃcult to extend those results to smooth domains
and the three dimensional cases.
This is a joint work with T. Suzuki and A. Mizutani, and a sequel of ourprevious work [4].

References[1] H. Brezis and W. Strauss: Semi-linear second-order elliptic equations in
L1 , J. Math. Soc. Japan, 25 (1973), 565–590.
[2] M. G. Crandall and T. Liggett, Generation of semi-groups of nonlinear
transformations on general Banach spaces, Amer. J. Math. 93 (1971),
265–293.
[3] R. Eymard, T. Gallouët and R. Herbin: Finite Volume Methods, Handbook Numer. Anal. VII, Elsevier, 2000, 713–1020.
[4] A. Mizutani, N. Saito and T. Suzuki: Finite element approximation fordegenerate parabolic equations. An application of nonlinear semigrouptheory, ESAIM: Math. Modell. Numer. Anal. 39 (2005), 755–780.

162

FEM analysis of the Stokes equations underboundary conditions of friction type
Takahito Kashiwabara
1) Graduate School of Mathematical Sciences, The University of Tokyo∗ e-mail: tkashiwa@ms.u-tokyo.ac.jp
Abstract1

Boundary conditions of friction type

We are concerned with the stationary Stokes equation in a two-dimensional
domain Ω ⊂ R2 , that is,−ν∆u + ∇p = f

in Ω,

div u = 0 in Ω,(1)

where ν, u, p, and f stand for a viscosity constant, velocity ﬁeld, pressure, andexternal force, respectively. As for the boundary, we assume that Γ := ∂Ω ispolygonal, that Γ is a union of two non-overlapping parts Γ0 , Γ1 , namely,Γ = Γ0 ∪ Γ1 ,

Γ0 ∩ Γ1 = ∅,and that Γ1 coincide with whole one side of the polygon Ω.
We associate with (1) the adhesive boundary condition u = 0 on Γ0 , whileon Γ1 we impose one of the following nonlinear boundary conditions:
un = 0,|στ | ≤ g,στ uτ + g|uτ | = 0,

(SBCF)

which is called the slip boundary condition of friction type, anduτ = 0,|σn | ≤ g,

σn un + g|un | = 0,

(LBCF)

which is called the leak boundary condition of friction type. Here, un (resp.
uτ ) is the normal (resp. tangential) component of u, στ (resp. σn ) is thetangential (resp. normal) components of the stress vector, and g is a given
positive number. In (SBCF) (resp. (LBCF)), if the magnitude of the tangential (resp. normal) stress is strictly less than g, then u = 0; if it reaches
the threshold g, then it can happen that uτ ̸= 0 (resp. un ̸= 0). This is why
the above boundary conditions are called of “friction type”.
In this study, ﬁnite element approximation of (1) under (SBCF) or (LBCF)
is examined. We propose approximate problems formulated by a variational
inequality, prove the unique existence of an approximate solution, establish
error analysis, and consider an Uzawa-type iteration method for numerical
implementation.
163

2Continuous and approximate problems

Due to space limitations, we restrict ourselves to (SBCF). Let {Th } be a
regular family of triangulations of the polygon Ω. Let Pk (T ) denote the setof the polynomial functions of degree ≤ k, deﬁned on T ∈ Th . We introduce:{}
Vn = v ∈ H 1 (Ω)2 v = 0 on Γ0 , vn = 0 on Γ1 ,◦{}Q = q ∈ L2 (Ω) (q, 1)0 = 0 ,}{Vnh = vh ∈ C(Ω)2 vh |T ∈ P2 (T )2 (∀T ∈ Th ), vh = 0 on Γ0 , vhn = 0 on Γ1 ,◦{}Qh = qh ∈ C(Ω) qh |T ∈ P1 (T ) (∀T ∈ Th ), (qh , 1)0 = 0 ,◦Now, the continuous and approximate problems, which seek (u, p) ∈ Vn × Q◦

and (uh , ph ) ∈ Vnh × Qh respectively, are formulated by the following variational inequalities:{a(u, v − u) + b(v − u, p) + j(v) − j(u) ≥ (f, v − u)0(∀v ∈ Vn ),(VI)◦b(u, q) = 0(∀q ∈Q).
{a(uh , vh −uh)+b(vh −uh , ph)+jh (vh)−jh (uh) ≥ (f, vh −uh)0 (∀vh∈Vnh ),(VIh )◦b(uh , qh ) = 0(∀qh ∈Qh ).
∫ ( ∂ui ∂uj )∑2
νHere, we deﬁne the bilinear forms a, b by a(u, v) = 2 i,j=1 Ω ∂xj + ∂xi
(
)
∫
∂vj
∂vi
dx, b(v, q) = − Ω divv q dx. The nonlinear term j, and its Simp+∂xj∂xi
son approximation jh , are deﬁned as∫j(v) = g|vτ | ds(v ∈ Vn ),jh (vh ) = g

Γ1
m
∑

)
|MiMi+1| (|vhτ (Mi )| + 4|vhτ (Mi+1/2 )| + |vhτ (Mi+1 )| (vh ∈ Vnh ),6i=1

where Mi ’s denote the vertices and the midpoints of two adjacent vertices,
contained in Γ1 and numbered in ascending order.
We will give details of our results, especially error analyses and numericalexamples, in the present talk. See also our preprint [1], where all of thetheoretical and numerical results are described in detail.

References[1] T. Kashiwabara: On a ﬁnite element approximation of the Stokes
problem under leak or slip boundary conditions of friction type,arXiv:1012.4982.
164

A mesh generator using a self-replicating
systemHirofumi NOTSU1),∗ , Masahiro YAMAGUCHI2) Daishin UEYAMA3)
1) Meiji Institute for Advanced Study of Mathematical Sciences
2) Graduate School of Science and Technology, Meiji University
3) School of Science and Technology, Meiji University∗ e-mail: notsu@meiji.ac.jp
Abstract
In this paper, a new triangular-mesh generator for numerical simulations,especially for ﬁnite element analyses, is presented.
For numerical simulations, meshes are necessary. However, it is not easy
to generate good meshes, because we need to add nodes, move the nodes
and change connections of triangles according to the shape of domain andrequired mesh size, where the details are not simple.
Now, complex patterns in nature have been modeled by nonlinear systems, which have been mathematically and computationally analyzed, e.g.,see [1]. We focus on so called self-replicating system. One of the simplest
self-replicating systems is the following Gray-Scott model [2], which is anonlinearly-coupled system of two partial diﬀerential equations and represents a chemical reaction.
∂u= Du ∆u − uv 2 + F (1 − u),
∂t
(1)
∂v

2

= Dv ∆v + uv + (F + k)v,∂t
where Du and Dv are diﬀusion constants and F and k are parameters. For
a parameter region, the system (1) give spot patterns. The properties of thepatterns to be created are as follows.
• The domain is ﬁlled with spots.
• The distances between spots are almost equal.
• The spots ﬁt the domain shape.
It is important that these properties are automatically realized in the system.
Using the system (1) and considering the spots as nodes, we generate a mesh.
Let Ω be a given domain and Γ be the boundary of Ω. For the given
domain Ω, what we need for computations is a mesh Th = {Kl }l , where Kl
is a triangle (2D) or a tetrahedra (3D) and h is a representative mesh size.
Let {Pi }i be nodes of the mesh Th . We propose a triangular-mesh generator
using the system (1), whose procedure is as follows.
1651. The domain Ω is given.
2. Artiﬁcial spots are put on the boundary Γ as an initial condition of thesystem (1).
3. The system (1) is numerically solved in an artiﬁcial rectangle domainΩ̃ (⊃ Ω) by a ﬁnite diﬀerence method (FDM).
4. A set of nodes {Pi }i is made by the created spots.
5. The Delaunay triangulation [3] is applied for the set of nodes and a setof triangles Th = {Kl }l is made.
We note that a rough computation is okay in the third process, and that themesh size h can be controlled by the ratio of diﬀusion constants Du and Dv .
Figure 1 shows snapshots of a mesh generation. It is observed that agood mesh has been generated for a complex domain. Since the procedure
is simple, and complex processes, e.g., smoothing, increasing and decreasing
nodes and so on, are solved by the self-replicating system, the mesh generator
is useful.

Figure 1: Snapshots of a mesh generation, a given domain (left), createdspots (center) and a generated mesh (right).

References[1] Y. Nishiura, Far-from-Equilibrium Dynamics, American MathematicalSociety, Providence, 2002.
[2] J. Pearson: Complex patterns in a simple system, Science, 261 (1993)
189–192.
[3] T. Taniguchi, Automatic mesh generation for FEM -Application of Delaunay triangulation, Morikita, Tokyo, 1992 (in Japanese).
166

A Remark on Numerical Instability ofComplex Inverse Laplace Transforms using
Multiple-Precision Arithmetic
Hiroshi FUJIWARA1),∗
1) Graduate School of Informatics, Kyoto University∗ e-mail: fujiwara@acs.i.kyoto-u.ac.jp
Abstract
We give a remark on numerical computation of the inverse Laplace transform based on the Bromwich integral. In particular numerical instability are discussed using multiple-precision arithmetic in Hosono’s inversion
method [1].
R∞
The Laplace transform Lf (s) = F (s) = 0 e−st f (t) dt appears as afundamental method in mechanics, engineering, mathematical science, andmathematical finance. And its inversion f = L−1 F also plays an importantrole. Under some assumptions, the Bromwich integral
Z c+iT1f (t) =
lim
etz F (z) dz,
2πi T →∞ c−iT
gives the inverse formula, and several numerical treatments have been proposed [2, 3]. Hosono has proposed a numerical scheme
!
µ
k−1
eσ X
1 X
f (t) ≈ fσ,N (t) =Aµ,ν Fk+ν ,Fn + k+1
t n=1
2
ν=0
where σ > 0 is an approximation parameter of the integral kernel etz , N
 =
1
1
n
(k, µ) is a truncation parameter,Fn = (−1)
Im F (zn ), zn = t {σ+ n − 2 πi},µ+1.
Aµ,µ = 1, and Aµ,ν−1 = Aµ,ν +ν
In the numerical treatments of the Bromwich integral, the numerical instability prevents reliable and accurate numerical computations. In Hosono’s
method the numerical solution is given by fσ,N = L−1
σ,N FN . The actual numerical computations are processed with floating-point arithmetic with precisionp, and the actual numerical result is
fσ,N,p = L−1
σ,N,p FN,p .
The difference |fσ,N − fσ,N,p | is caused by floating-point arithmetic, and it iscalled a rounding error in the precision p.
167

We propose the use of multiple-precision arithmetic [4] for quantitative
analysis of rounding errors. Figure 1 shows quantity of rounding errors by
fσ,N,p (t) − fσ,N (t),ulp fσ,N,p (t)rounding error [ULP] (logarithmic scale)where ulp(x) is the unit in the last place of a floating-point value x. From thefigure, the magnitude of rounding errors in the numerical solution is approximately 1016 for σ = 40, and the standard double precision arithmetic, whichhas 16 decimal digits precision, is not enough for the reliable computations.
1e+25

σ = 40σ = 20σ = 10

1e+20

1e+15

1e+10

100000

1

0

5

10

15

20
t

25

30

35

40

Figure 1: Rounding Errors in Numerical Solution by Hosono’s Method
The truncation error |f − fσ,N | tends to zero as σ tends to infinity in
Hosono’s method. We need large σ to reduce truncation errors, which leads
the numerical instability. Multiple-precision arithmetic is effective for accurate and reliable numerical computations.

References[1] T. Hosono: Numerical Inversion of Laplace Transform (in Japanese),
IEEJ Trans. FM, 99 (1979), 494–50.
[2] A. M. Cohen: Numerical methods for Laplace transform inversion,
Springer, 2007.
[3] D. Sheen, I. H. Sloan, and V. Thomée: A parallel method for time discretization of parabolic equations based on Laplace transformation andquadrature, IMA J. Num. Anal. 23 (2002), 269–299.
[4] exflib: http://www-an.acs.i.kyoto-u.ac.jp/~ fujiwara/exflib
168

Numerical Laplace inversion using multi-precision
Jiwoon Kim∗

Dongwoo Sheen†

Abstract
Numerical inversion of Laplace transform has attracted intensive attention fromscientists and engineers as well as mathematicians. Various numerical methods havebeen devised. In particular, we review several efficient methods based on deformed
contour that have been developed for about ten years. A comparative study using
multi-precision will be exposed. We examine the effects of roundoff errors, computation
is carried out both in double precision and a multi-precision, the latter which provides
better understanding of the numerical Laplace inversion algorithms. In this talk wealso review and compare the most efficient methods for numerical Laplace inversion,in particular, using multi-precision. In this study, several methods have been classified
and interpreted as the contents related to three issues: (i) the choice of contour, (ii)its parameterization, and (iii) numerical quadrature rule.

∗

Department of Mathematics, Seoul National University, Seoul 151-747, Korea; E-mail: jiwoonkim83@gmail.com†
the same address; also, Interdisciplinary Program in Computational Science & Technology, Seoul National University, Seoul 151–747, Korea; E-mail: sheen@snu.ac.kr

1169CPU-GPU Hybrid Approaches inMultifrontal Methods for Large and Sparse
Linear System
Chenhan D. Yu1),∗ , Weichung Wang1) , Dan’l Pierce2)
1) Department of Mathematics, National Taiwan University, 10617, Taiwan2) Access Analytics Intl, Redmond, Washington 98073, USA∗ E-mail: b94201001@gmail.com
Abstract
Solving large-scale sparse linear systems Ax = b is at the heart of various scientiﬁc and engineering computations. To solve these problems, direct
methods relying on symbolic and numerical factorizations have been used for
decades and remain popular choices, as evidenced by the many well-designed
algorithms and well-developed software packages that are available. Most
recently, a growing number of computing programs have tried to use graphic
processing units (GPU) to accelerate the computations, due to its highly
parallel computational capability. We intend to investigate how we can integrate a GPU with a central processing unit (CPU) system to accelerate thedirect solving process for sparse linear systems.
Among various direct methods, we focus on the multifrontal method particularly [1, 2]. The choice of multifrontal method is mainly motivated by
the following observation. A multifrontal method uses a factorization tree totransform a large sparse linear system problem into many operations involving smaller dense matrices. These dense operations are promising candidates
that can be performed on a GPU to gain possible time savings without downgrading the accuracy. As all of the frontal matrices are dense and the main
computational tasks involving these dense frontal matrices are actually level
2 and level3 Basic Linear Algebra Subprograms (BLAS) computations, we
can accelerate the whole multifrontal process if we can accelerate the BLAS
computations via a GPU.
A GPU contains a large number of computational cores on a single board.
These cores can work simultaneously and thus justify the ideas to accelerate
the BLAS and other operations in a multifrontal method via a GPU. However, it is important to realize that a GPU is not a standalone platform, but
a device driven by a CPU. Consequently, the data must be transferred between the host CPU and the device GPU. The communication cost between
the CPU and GPU thus plays a critical role in the overall performance.
For unsymmetric linear systems, we analyze the multifrontal method fromboth an algorithmic and implementation perspective to see how a GPU, in170particular the NVIDIA Fermi series, can be used to accelerate the computations. Our main accelerating strategies include (i) performing BLAS on
both CPU and GPU, (ii) improving the communication eﬃciency between
the CPU and GPU by using page-locked memory and asynchronous memory
copy, and (iii) a modiﬁed algorithm that reuses the memory between diﬀerent GPU tasks and sets thresholds to determine whether a certain tasks be
performed on GPU.
For symmetric positive deﬁnite systems, the step-by-step operations of
the numerical factorization can be exactly predicted after the symbolic factorization is done. However, data communication between CPU and GPU
remains a critical issue that can aﬀect the overall performance signiﬁcantly.
We analyze the symmetric multifrontal method to ﬁnd the minimal communication requirements between CPU and GPU for symbolic factorization,coeﬃcient matrix, numerical factorization and triangular solver. Besides,
we propose a timing-optimized algorithm that further improve the overall
performance. In the algorithm, some data communications are hidden byoverlapping the data transfers with other operations. In addition, we show
how shorter overall computational time can be achieved by setting thresholds to let small frontals be factored in CPU, even extra data communication
traﬃcs are introduced.
The proposed acceleration strategies are implemented by modifying UMFPACK (http : //www.cise.uf l.edu/research/sparse/umf pack/) and TAUCS(http : //www.tau.ac.il/ stoledo/taucs/) which is an unsymmetric and symmetric multifrontal linear system solver, respectively. Numerical results fora series of three-dimensional Poisson problems and standard large test problems will be presented. The results show that the CPU-GPU hybrid approach
can accelerate the unsymmetric and symmetric multifrontal solver, especiallyfor computationally expensive problems.

References[1] T. Davis, A column pre-ordering strategy for the unsymmetric-pattern
multifrontal method, ACM Transactions on Mathematical Software(TOMS) 30 (2) (2004) 195.
[2] J.Liu,The multifrontal method for sparse matrix solution: Theory andpractice, SiamReview 34 (1) (1992) 82109.

171

Design and implementation of amultiple-precision system on GPU
Shugo Manabe1),∗1) Graduate School of Informatics, Kyoto University∗ e-mail: shugo@acs.i.kyoto-u.ac.jp
Abstract
In the present research, we design and implement a multiple-precisionsystem which is useful for parallel computation.
In scientific computation, there are some hard problems to obtain desired
accuracy of the results due to accumulation of rounding errors. We propose
adopting multiple-precision system, but arithmetic speed is slow on such
systems. In scientific computation, vector operations usually occupy mostof the computation time. GPU (Graphics Programming Unit) has highly
parallel structure. In the current research, we design and implement themultiple-precision system on GPU.
GPU is a specialized microprocessor that offloads and accelerates 3D or
2D grahics rendering from CPU (Central Processing Unit). Modern GPU is
not only a powerful graphics engine but also a highly parallel programmableprocessor. Nowadays application of GPU to scientific computation develops.
GPU has a large number of cores. For example, Intel Core i7-950 (CPU) has
only 4 cores, on the other hand Geforce GTX 480 (GPU) has 480 cores.
Now we have implemented 1) vector operations, 2) decimal input andoutput, 3) type conversion, and 4) summation and inner product on the proposed system. In particular, 1) and 4) are optimized for GPU. The following
are typical vector operations.
Z[i] = X[i] + Y [i](i = 0, 1, . . . , N − 1)
Z[i] = X[i] − Y [i](i = 0, 1, . . . , N − 1)
Z[i] = X[i] × Y [i](i = 0, 1, . . . , N − 1)
Z[i] = X[i] ÷ Y [i](i = 0, 1, . . . , N − 1)
where X[i], Y [i], and Z[i] are multiple-precision numbers.

(1)
(2)
(3)
(4)

We do not have to care how to use GPU. For example, if we compute (1),
we only write
vectZ = vectX + vectY;172

on the C++ programming language. We can use the proposed system without knowledge of GPU.
Figure 1 and Figure 2 compare the proposed system with exflib (multipleprecision arithmetic library)[1, 2]. On the proposed system, as the numberof elements of vector gets larger, you can obtain more performance.
×106 operations per sec

proposed system

exflibN=106
5
Geforce GTX 480 N=104
N=10
N=103
Athlon64 2.20GHz
Core i7 2.66GHz

3833.5(×130.3)
3332.8(×113.3)
1875.3(×63.8)
218.3(×7.4)
29.4(×1.00)
14.8(×0.50)- FastFigure 1: comparison of speed of addition with 100 decimal precision

×106 operations per sec

proposed system

exflibN=106
5
Geforce GTX 480 N=104
N=10
N=103

307.7(×40.9)
303.7(×40.4)
250.1(×33.3)
54.4(×16.9)

Athlon64 2.20GHz
Core i7 2.66GHz

7.52 (×1.00)
3.22 (×0.50)- FastFigure 2: comparison of speed of multiplication with 100 decimal precision

References[1] exflib home page; http://www-an.acs.i.kyoto-u.ac.jp/∼ fujiwara/
exflib/index.html[2] Hiroshi Fujiwara : Design and implementation of a fast multipleprecision environment for scientific computation (in Japanese), Jourmnal of the Society of Instrument and Control Engineers , vol. 49, No. 5(2010) pp. 285–290.
173

On Optimality Theorems for Robust
Multiobjective Optimization Problems
Gue Myung Lee
Department of Applied Mathematics, Pukyong National University,
608-737, Korea
e-mail: gmlee@pknu.ac.kr
Abstract
Recently, many authors have studied roust (scalar) optimization problems([1-12]). In this talk, the robust approach (the worst-case approach) for the
standard multibjective optimization problem (MP) is considered. We define
three kinds of robust efficient solutions for the robust counterpart (RMP)of the uncertain multiobjective optimization problem (UMP) which consistsof more than two objective functions with uncertainty data and constraintfunctions with uncertainty data. We give necessary optimality theorems forweakly and properly robust efficient solutions for (RMP).
This talk is based upon the paper ”Robust Multiobjective Optimization”
written by Daishi Kuroiwa and Gue Myung Lee.

References[1] A. Beck and A. Ben-Tal, Duality in robust optimization: primal worst
equals dual best, Oper. Res. Lett., 37 (2009), 1-6.
[2] A. Ben-Tal and A. Nemirovski, Robust optimization—methodology andapplications, Math. Program., Ser B, 92 (2002), 453-480.
[3] A. Ben-Tal, L.E. Ghaoui and A. Nemirovski, Robust Optimization,
Princeton Series in Applied Mathematics, 2009.
[4] A. Ben-Tal and A. Nemirovski, A selected topics in robust convex optimization, Math. Progr. B, 112 (2008), 125-158.
[5] D. Bertsimas and D. Brown, Constructing uncertainty sets for robust
linear optimization, Oper. Res., 57 (2009), 1483-1495.
[6] D. Bertsimas, D. Pachamanova and M. Sim, Robust linear optimization
under general norms, Oper. Res. Lett., 32 (2004), 510-516.

174

[7] V. Jeyakumar and G. Li, Characterizing robust set containments andsolutions of uncertain linear programs without qualifications, Oper. Res.
Lett., 38(2010), 188-194.
[8] V. Jeyakumar and G. Y. Li, Robust Farkas’ lemma for uncertain linear
systems with applications, Positivity, DOI 10. 1007/s11117-010-0078-4.
[9] V. Jeyakumar and G. Y. Li, Strong duality in robust convex programming: complete characterizations, SIAM J. Optim., 20(2010), 33843407.
[10] V. Jeyakumar, G. Y. Li and G. M. Lee, A robust von Neumann minimax
theorem for zero-sum games under bounded payoff uncertainty, Oper.
Res. Lett., 39(2011), 109-114.
[11] V. Jeyakumar, G. Y. Li and G. M. Lee, Robust duality for generalized convex programming problems under data uncertainty, to appear
in Nonlinear Anal..
[12] G. Y. Li, V. Jeyakumar and G. M. Lee, Robust conjugate duality forconvex optimization under uncertainty with application to data classification, Nonlinear Anal., 74(2011), 2327-2341.

175

The Resolvent average on symmetric cones
Sangho Kum1),∗ , Yongdo Lim2)
1) Department of Mathematics Education, Chungbuk National University,2) Department of Mathematics, Kyungpook National University∗ e-mail: shkum@cbnu.ac.kr
Abstract
Recently Bauschke et al. introduced a very interesting and new notion ofproximal average in the context of convex analysis, and studied this subject
systemically in [3–7] from various viewpoints. In addition, this new concept
was applied to positive semidefinite matrices under the name of resolvent
average, and basic properties of the resolvent average are successfully established by themselves from a totally different view and techniques of convex
analysis rather than the classical matrix analysis [8]. Inspired by their works
and the well-known fact that the convex cone of positive definite matrices isa typical example of a symmetric cone (self-dual homogeneous convex cone),we study the resolvent average on symmetric cones, and derive corresponding
results in a different manner based on a purely Jordan algebraic technique.
By the way, the functional f (x) = − ln det(x) on a symmetric cone isa self-scaled barrier. Indeed, this logarithmic barrier function f (x) lies at
the heart of interior point methods in optimization over symmetric cones[1, 12, 13, 14]. In regard to this fact, we provide a variational characterizationof the resolvent average as a unique minimizer of a Bregman distance inducedby the parameterized logarithmic barrier functional fµ (x) = −ln det(x + µe).

References[1] M. Baes, Convexity and differentiability properties of spectral functions
and spectral mappings on Euclidean Jordan algebras, Linear AlgebraAppl. 422 (2007), 664–700.
[2] H. H. Bauschke, E. Matouskova and S. Reich, Projection and proximal point methods: convergence results and counterexamples, Nonlinear
Anal. TMA 56 (2004), 715–738.
[3] H. H. Bauschke, Y. Lucet and X. Wang, Primal-dual symmetric intrinsic
methods for finding antiderivatives of cyclically monotone operators,
SIAM J. Control Optim. 46 (2007), 2031–2051.
176

[4] H. H. Bauschke, Y. Lucet, and M. Trienis, How to transform one convex
function continuously into another, SIAM Rev. 50 (2008), 115–132.
[5] H. H. Bauschke, R. Goebel, Y. Lucet, and X. Wang, The proximal
average: basic theory, SIAM J. Optim. 19 (2008), 766–785.
[6] H. H. Bauschke and X. Wang, The kernel average for two convex functions and its applications to the extension and representation of monotone operators, Trans. Amer. Math. Soc. 361 (2009), 5947–5965.
[7] H. H. Bauschke, Y. Lucet and M. Trienis, The piecewise linear-quadraticmodel for computational convex analysis, Comput. Optim. Appl. 43(2009), 95–118.
[8] H. H. Bauschke, S. M. Moffat and X. Wang, The Resolvent average forpositive semidefinite matrices, Linear Algebra Appl. 432 (2010), 1757–
1771.
[9] S. Kim, J. Lawson and Y. Lim, The matrix geometric mean of parameterized, weighted arithmetic and harmonic means, submitted.
[10] J. Faraut and A. Koranyi, Analysis on Symmetric Cones, Clarendon
Press, Oxford, 1994.
[11] Y. Lim, Birkhoff formula for conformal compressions of symmetric cones,
Amer. J. Math. 125 (2003), 167-182.
[12] Y. E. Nesterov and M. J. Todd, Self-scaled barriers and interior-point
methods for convex programming, Math. Oper. Res. 22 (1997), 1-42.
[13] Y. E. Nesterov and M. J. Todd, On the Riemannian geometry definedby self-concordant barriers and interior-point methods, Found. Comput.
Math. 2 (2002), 333-361.
[14] Y. E. Nesterov and A. Nemirovski, Primal central paths and Riemannian
distances for convex sets, Found. Comput. Math. 8 (2008), 533–560.

177

Optimizing Latin Hypercube Designs by
Particle Swarm with GPU Acceleration
Dai-Ni Hsieh1),∗ , Ray-Bing Chen2) , Ying Hung3) , Weichung Wang1)
1) Department of Mathematics, National Taiwan University, Taiwan2) Department of Statistics, National Cheng Kung University, Taiwan
3) Department of Statistics and Biostatistics, Rutgers University, USA∗ E-mail: r97221050@ntu.edu.tw
Abstract
Computer and physical experiments of many kinds are important tools
in scientific researches. As many of these experiments are time consuming,
it is important to design the experiment carefully with small number of experimental points. In order to explore an experimental region efficiently, theexperimental design we used in a computer experiment is required to satisfy
certain properties like space-filling and non-collapsing. Based on these two
properties, a space-filling Latin hypercube design (LHD) is an appropriate
and popular choice.
Formally, an n-run and k-factor LHD is an n by k matrix. Each row of
the matrix corresponds to a design point, and each column is a permutationof {1, 2, . . . , n}. Note that the total number of n-run and k-factor LHDs is(n!)k−1 , which grows fast as n or k increase. Although a LHD is guaranteed to
be non-collapsing, a randomly generated LHD is usually not space-filling. To
satisfy the space-filling property, Morris and Mitchell [1] suggested applying
the maximin criterion to LHDs, i.e. maximizes the minimal pairwise distanceof design points. Furthermore, they proposed the extended maximin criterion
to break ties among optimal maximin LHDs . In order to use the extended
maximin criterion in a easier way, they introduced the φp criterion,


n−1
X

φp = 

n
X

1

p

 ,
d−p
ij

(1)

i=1 j=i+1

where p is a positive integer and dij is the distance between the i-th and j-th
design points. When p is large enough, φp ranks LHDs in the same way as
the extended maximin criterion does. In this paper, we search the optimal
LHD by minimizing φp with respect to all possible LHDs.
The main challenge of this work arises from the huge number of feasible
points in the discrete optimization problem min φp . To tackle the problem,
we modify the particle swarm optimization (PSO) [2] framework to fit itinto the LHD optimization problem. PSO is a population based stochastic
178

heuristic methodology and is capable of solving high-dimensional optimization problems with multiple optima. The key idea of PSO is that all particles
cooperate to find the optimum in the search domain through social interaction. However, the original PSO is to solve continuous optimization problems.
Following the structure of PSO, we alter the step of position update such thatthe movement of a particle is from one LHD to another LHD. Moreover, we
use graphic processing unit (GPU) to accelerate the proposed PSO.
To assess the quality of LHDs we found, the squared minimum distance of optimal LHDs generated by PSO and other approaches are presented in Table 1. In the table, the column UB shows the upper bounds
obtained by van Dam et al. [3]. For the column Rand, we show the bestvalues of randomly generated LHDs whose total number is 100 times thenumber of LHDs visited in PSO. The column Web are the results fromhttp://www.spacefillingdesigns.nl, which contains optimal maximin LHDsof several papers. Comparing with random LHDs, PSO always gets better
results. If we compare with the results from the website, PSO is capable ofachieving the same value and even making some improvements. For others
that are not available from the website, the results obtained by PSO are nearto the upper bounds. In addition to the better squared minimum distance,
we usually get more than 20X speedup on a NVIDIA Telsa C2070 GPU
comparing to the CPU version that is run on a HP BL460cG6 workstationequipped with Intel Xeon X5570 2.93GHz CPU. These results together withthe acceleration of GPU suggest that the proposed PSO is promising to solve
even larger problems.
Table 1: Optimal results comparison. (A higher value is better.)
n

k

UB

Rand

Web

PSO

n

k

UB

Rand

Web

PSO

6
6
6
6
6
6

10
11
12
13
14
15

70
77
84
91
98
105

65
72
79
85
92
99

68
74
82
89
N/A
N/A

68
75
84
89
96
104

7
7
7
7
7
7

10
11
12
13
14
15

93
102
112
121
130
140

83
91
101
109
118
128

90
N/A
N/A
N/A
N/A
N/A

90
100
109
119
128
138

References[1] M.D. Morris, T.J. Mitchell: Exploratory designs for computational experiments, J. Statist. Plann. Inference, 43 (1995), 381–402.
[2] J. Kennedy, R.C. Eberhart: Particle swarm optimization, Proc. IEEE
international conference on neural networks, 4 (1995), 1942–1948.
[3] E.R. van Dam, G. Rennen, B. Husslage: Bounds for maximin Latin
hypercube designs, Oper. Res., 57 (2009), 595–608.
179

Maximum Entropy Approach For Optimal Repairable Mx/G/1
Queue With Bernoulli Feedback And Setup
M. JainIndian Institute of Technology, Roorkee-247667 (India)       
Abstract
Queueing theory was developed to predict the behavior of various service systems
wherein the server may experience unpredictable breakdowns. Therefore queueing
systems with unreliable servers reflect more practical situations that we face in real timesystems. Bernoulli feedback queueing models have been developed to analyze many
stochastic service systems such as telecommunication systems, manufacturing andproduction systems, packet transmission, etc. In such type of systems, the customers are
fed back again to the head of the queue in case if he is not satisfied with the service
provided by the server in the first attempt. This process is continued until the customer isfully satisfied and leaves the system after the completion of the service.
The model developed in this paper has been designed to analyze the behavior of
Mx/G/1 queueing system with balking. The server initiates the service only if there are N
customers accumulated in the system. To make the model more versatile, the assumptionsof unreliable server, general phase repair, Bernoulli feedback, startup and setup times aretaken into consideration. Using the supplementary variable technique and the generating
function method, the stationary distribution of the queue length and various otherqueueing characteristics are established.
Using principle of maximum entropy (PME), one can propose a feasible method
for approximately analyzing many complex queueing problems. The idea behind thisapproach is to develop a maximum entropy analysis which in turn provides the widest
probability distribution subject to the known mean constraints. The PME used for
estimating the probabilistic information measures has been employed to obtain queue
size distribution of various queueing systems in different frameworks. Using the
maximum entropy approach, we intend to determine the approximate results for thesteady-state probability distributions of the queue length. The method of Lagrange’s
undetermined multipliers to maximize the entropy function subject to the relevant
constraints is used. Moreover, we derive the approximate formulae for the expectedwaiting time of an arbitrary customer in the queue. A comparative study between theexact results obtained using supplementary variable technique and the approximate
results obtained using maximum entropy approach has been facilitated by taking
numerical illustration. The comparative analysis facilitate an insight to the practitioners
and system engineers to make appropriate choice of the service time distribution andrepair time distribution, which give the minimum expected waiting time. The numericalillustration provided indicate that the PME is accurate enough for practical purposes andprovides a useful method to cope up with many complex queueing congestion situations
which are not easily explored by using classical analytical techniques.

180

Application of the Bees Algorithm to Multi-Objective
Optimization Engineering Problems
A.Ghanbarzadeh1),* , A.Moradi2) , A.Mirzakhani Nafchi3)
1) School of Mechanical Engineering, Shahid Chamran University, Ahvaz,

IRAN2) Graduate School of Mechanical Engineering, Shahid Chamran University,
Ahvaz, IRAN
3) Graduate School of Mechanical Engineering, Shahid Chamran University,
Ahvaz, IRAN
*Ghanbarzadeh.A@scu.ac.ir, A.Moradi64@gmail.com,
Amin.Mirzakhani@yahoo.com
Abstract
This paper describes application of the Bees Algorithm [1] to multi-objective
optimization problems. The Bees Algorithm is a search procedure inspired by theforaging behavior of honey bees [2]. In the present paper, two standard mechanicaldesign problems, the design of a gear train design and the design of coil springs,
were used to benchmark the Bees Algorithm to compare with other optimization
techniques.
Multi-objective optimization using the Bees Algorithm is developed in order toobtain a set of geometric design parameters, lead to find optimum solution. Multiobjective optimization procedure yields a set of non-dominated solutions, called
Pareto optimal set, each of which is a trade-off between objectives and can beselected by the user, regarding the application and the project’s limits. The
presented work takes care of numerous geometric parameters in the presence of
logical constraints. Study results showed that the approach using the Bees
Algorithm, found better solutions than traditional calculus-based algorithm, BFGS.
The coil spring design problem consists of minimizing the weight of atension/compression spring subject to constraints on shear stress, surge frequency
and minimum deflection. The design variables are the mean coil diameter D ; thewire diameter d, and the number of active coils N [3].
The objectives of the gear train design is to find the number of teeth in each of
the four gears so as to minimize (i) the error between the obtained gear ratio and arequired gear ratio of 1/6.931 and (ii) the maximum size of any of the four gears.
Since the number of teeth must be integers, all four variables are strictly integers
[3].

181

The Bees Algorithm has been successfully applied to various benchmark andstandard engineering optimization problems. Numerical results reveal that the
proposed algorithm can find better solutions when compared to other methods andis a powerful search algorithm for various engineering optimization problems.

References[1] Pham D.T., Ghanbarzadeh A., Koc E., Otri S., Rahim, S., and Zaidi M., 2005.
“Technical Note: The Bees Algorithm”. Technical Report No MEC 0501,
Manufacturing Engineering Centre, Cardiff University: Cardiff.
[2] Von Frisch, K. Bees: their vision, chemical senses and language, revised
edition, 1976 (Cornell University Press,Ithaca, NY).
[3] Deb, K., Pratap, A., and Moitra, S., 2000. “Mechanical Component Designfor Multiple Objectives Using Elitist Non-Dominated Sorting GA”. Technical
Report No 200002, Indian Institute of Technology: Kanpur, India. pp.10.

182

A quadrilateral Morley element forbiharmonic equations
Chunjae Park1),∗ , Dongwoo Sheen2)
1) Department of Mathematics, Konkuk University2) Department of Mathematics, Seoul National University∗ e-mail: cjpark@konkuk.ac.kr
Abstract
We consider the biharmonic problem in a polygonal domain Ω ⊂ R2 ,∆2 u = f
u = 0
∂u
= 0
∂ν

in Ω,
on ∂Ω,on ∂Ω.

(1)

2
The variational formulation of (1) is to find a solution u ∈ H00(Ω) such thata(u, v) = (f, v),

2∀v ∈ H00
(Ω),

(2)

where

∂u
= 0 on ∂Ω},
∂ν22and the bilinear form a : H00(Ω) × H00(Ω) → R which is defined byZ
2a(u, v) =uxx vxx + 2uxy vxy + uyy vyy dxdy, ∀u, v ∈ H00
(Ω).
2
H00
(Ω) = {u ∈ H 2 (Ω) | u =Ω2(Ω) such as ArThere are several conforming finite element spaces in H00
gyris, Bell, Hsieh-Clough-Tocher, Singular Zienkiewicz triangles and Fraeijs
2de Veubeke-Sander quadrilateral [2]. To be in H00(Ω) is to be of class C 1
for a finite element space in practice. In order to solve (2) with a finite
element space of lower regularity, some mixed methods and nonconforming
finite element methods have been developed [1].
For given triangulation T of Ω, define a space W (T ) asZ∂v
2
W (T ) = {v ∈ L (Ω) v(P ),dσ are continuousE ∂ν
at all vertices P and edges E in T }.

The Morley finite element space is one of the simplest nonconforming finite
element spaces with the simplicial triangulation T S of Ω. It is a subspace183

of W (T S ) and its element is piecewise quadratic. Some rectangular Morley finite element spaces, that is, subspaces of W (T R ) to solve (2) with arectangular triangulation T R have been designed [3].
Let T Q be a triangulation of Ω into convex quadrilaterals. In this talk,we will introduce a quadrilateral Morley finite element space in W (T Q ). For
each quadrilateral Q ∈ T Q , our finite element space belongs toM
P 2 (Q)< L1 L3 L13 , L2 L4 L24 >,where Li is a linear function which vanishes at edge Ei of Q for i = 1, · · · , 4
and Li i+2 is also linear which vanishes at the bimedian between Ei , Ei+2 for
i = 1, 2. Ei , i = 1, · · · , 4 are numbered counterclockwisely as in Figure 1.
L2 =0
E2

E1
E3

L3=0

E4

L13 =0

L24=0

L4=0
L1=0

Figure 1: Name convention in Q
We will prove the values at vertices and integrals of normal derivatives
over edges are the degrees of freedom of the proposed element. An optimalorder of convergence is analyzed and several numerical tests confirm it.

References[1] D. N. Arnold, F. Brezzi : Mixed and nonconforming finite element
methods: implementation, postprocessing and error estimates, RAIRO
Modél. Math. Anal. Numér., 19 (1985), 7–32.
[2] P. G. Ciarlet, J. L. Lions: Handbook of Numerical Analysis, Vol. II,
North Holland, 1991.
[3] M. Wang, Z. Shi, J. Xu : Some n-Rectangle Nonconforming Elements
for Fourth Order Elliptic Equations, Journal of Computational Mathematics, 25 (2007), 408–420.
184

Applying Snapback Repellers in Ecology
Shu-Ming Chang
Department of Applied Mathematics, National Chiao Tung University
e-mail: smchang@math.nctu.edu.tw
Abstract
In ecology, Satake’s generalized resource budget model [6] that modiﬁed
from Isagi’s resource budget model [2], Satake and Iwasa illustrated by computing the positive Lyapunov exponent that if the depletion coeﬃcient isgreater than one, then the system is chaotic. However, a positive Lyapunov
exponent implies only sensitivity in Devaney’s chaos [1, 5]. Therefore, thiswork presents mathematical viewpoints and numerical analysis on Satake’s
generalized resource budget model, to rigorously prove that the generalized
resource budget model is chaotic in Devaney’s sense by using the snapback
repeller theory [4, 7] and the topological entropy theory [3]. Moreover, thiswork also investigates that the behaviors are diﬀerent between positive odd
depletion coeﬃcients and positive even depletion coeﬃcients under numericalcomputations.

References[1] R. L. Devaney. An Introduction to Chaotic Dynamical Systems, Second
Edition. Addison-Wesley, Redwood City, Canada, 1989.
[2] Y. Isagi, K. Sugimura, A. Sumida, and H. Ito. How does masting happen
and synchronize? J. theor. Biol., 187:231–239, 1997.
[3] D. Kwietniak and M. Misiurewicz. Exact devaney chaos and entropy.
Qualitative Theory Dynamical Systems, 6(95):169–179, 2005.
[4] F. R. Marotto. On redeﬁning a snap-back repeller. Chaos, Solitons andFractals, 25:25–28, 2005.
[5] C. Robinson. Dynamical Systems: Stability, Symbolic Dynamics, andChaos, Second Edition. CRC, Boca Raton, Florida, 1998.
[6] A. Satake and Y. Iwasa. Pollen coupling of forest trees: Forming synchronized and periodic reproduction out of chaos. J. theor. Biol., 203:63–
84, 2000.
[7] Y. Shi and P. Yu. Chaos induced by regular snap-back repellers. J.
Math. Anal. Appl., 337:1480–1494, 2008.
185

Verified solutions of linear systems on GPU
Yusuke Morikura1),∗ , Katsuhisa Ozaki2) , Shin’ichi Oishi3)
1) Graduate School of Fundamental Science and Engineering,
Waseda University2) Department of Mathematical Sciences,
Shibaura Institute of Technology & JST, CREST
3) Faculty of Science and Engineering, Waseda University & JST, CREST∗ e-mail: m.myusuke@suou.waseda.jp
Abstract
This talk is concerned with accuracy of numerical solutions of liner systems, where numerical solutions mean computed results by ﬂoating-point
arithmetic deﬁned by the IEEE 754 standard [1]. Due to accumulation ofrounding errors, the numerical solution may be inaccurate. The aim of our research is to obtain an error bound for a numerical solution of a linear system.
This topic is frequently discussed in the ﬁeld of veriﬁed numerical computations. Let x̃ be an approximate solution of a linear system Ax = b where
A is a real n × n matrix and b is a real n-vector. Let R be an approximate
inverse matrix of A and I be the n × n identity matrix. If kRA − Ik∞ < 1,then the matrix A is nonsingular andkx̃ − A−1 bk∞ ≤kR(Ax̃ − b)k∞.
1 − kRA − Ik∞(1)

Based on (1), several algorithms obtaining an upper bound of the error havebeen developed [2, 3, 4].
Nowadays, GPU (Graphics Processing Unit) is used for not only acceleration of building of images but also for numerical computations since theperformance is very high. For example, Fermi from Nvidia Corporation has
448 cores and peak performance is 515 Gﬂops for double precision arithmetic. Basically, it was diﬃcult for a beginner for high performance computing to exploit GPU eﬃciently. Recently, useful toolboxes for GPU has beendeveloped, for example, CUDA, JACKET and so forth. Therefore, GPU
will be widely used. However, there are two problems of implementation ofalgorithms for veriﬁed numerical computations on GPU. Rounding modesdeﬁned by the IEEE 754 standard (rounding to nearest, rounding upward,
rounding downward) are useful for veriﬁed numerical computations. For example, they are used for interval arithmetic. If we use routines supported
in such libraries, then the rounding mode cannot be switched dynamically.
186

An algorithm on GPU must output the error bounds without switches of therounding modes, namely, only rounding to nearest is used. Next, amount ofworking memory is little comparing to CPU’s space, for example, C2050 and
C2070 by Nvidia Corporation have 3 GBytes and 6 Gbytes working memory,
respectively. Therefore, computable size of the matrix is limited.
We develop an algorithm for the veriﬁcation of linear systems which is
specialized for GPU. Our algorithm only requires ﬂoating-point computations
with the rounding to nearest mode. Moreover, applying block algorithms, weovercome these problems. Finally numerical result will be shown to illustrate
the eﬃciency of the proposed algorithm.

References[1] ANSI/IEEE: IEEE Standard for Binary Floating-Point Arithmetic:
ANSI/IEEE Std 754-1985, New York, IEEE, 1985 .
[2] S. Oishi, S. M. Rump: Fast Veriﬁcation Solutions of Matrix Eqiations,
Numer.Math.,90:4(2002),755-773.
[3] T. Ogita, S. M. Rump, S. Oishi: Veriﬁed Solution of Liners Systems
without Directed Rounding,11th GAMM-IMACS International Symposium on Scientiﬁc Computing. Computer Arithmetic, and Validated Numerics, Fukuoka,Japan,2004[4] K. Ozaki, T. Ogita, S. Miyajima, S. Oishi, S. M. Rump: A Method of
Obtaining Veriﬁed Solutions for Linear Systems suited for Java, Journalof Computational and Applied Mathematics, 199:2(2007), pp. 337-344.

187

Speaker Index
Name

Page no.

Session id

Afshin Ghanbarzadeh

21 181[29-4-C]

Akitoshi Takayasu

12 33[27-3-A]

Akiyasu Tomoeda

15 78[28-3-C]

Alex Chang

17 105[28-4-C]

Andrea Bertozzi

18 121[29-1-A]

Buyung Sarita

16 91[28-3-D]

Carolina Mendoza

15 80[28-3-C]

Chang Yong Han

19 137[29-2-D]

Changtong Luo

13 60[27-4-D]

Chenhan Yu

21 170[29-4-B]

Chien-Hong Cho

13 49[27-4-C]

Chin-Tien Wu

18 124[29-2-B]

Ching-An Wang

13 44[27-4-C]

Chiun-Chang Lee

12 32[27-3-A]

Chunjae Park

22 183[29-4-D]

Dai-Ni Hsieh

21 178[29-4-C]David Ni

18 122[29-2-B]

Dongwoo Sheen

21 169[29-4-B]

Edi Cahyono

16 90[28-3-D]

Elliott Ginder

16 96[28-4-B]

Eun-Jae Park

14 66

[28-2-A(1)]G Sharma

17 112[28-4-C]

Guanyu Zhou

20 159[29-3-D]

Gue Myung Lee

21 174[29-4-C]( to be continued )188

Speaker Index
Name

Page no.

Session id

HUI JI

11 27[27-2-A(2)]Haiwei Sun

16 101[28-4-B]

Hendra Gunawan

13 51[27-4-C]

Hidenori Yasuda

15 70[28-3-B]

Hirofumi Notsu

20 165[29-3-D]

Hiroshi Fujiwara

21 167[29-4-B]

Hiroshi Kokubu

16 85[28-3-D]

Hiroshi Suito

17 103[28-4-C]

Hisashi Okamoto

15 84[28-3-C]

Ikha Magdalena

17 110[28-4-C]

Inseok Yang

17 119[28-4-D]

Kaname Matsue

16 86[28-3-D]

Katsuhisa Ozaki

12 42[27-4-B]

Keiji Miura

17 115[28-4-D]

Keita Iida

15 82

[28-3-C]

Kenta Kobayashi

18 131[29-2-C]

Kil H Kwon

18 128[29-2-B]

Kuo-Chih Hung19 145[29-3-B]

Kwok Kin Wong

13 53

[27-4-C]

Lei Du12 30[27-3-A]

Leo Hari Wiryanto

13 55[27-4-D]

Madhu Jain

21 180[29-4-C]

Masataka Kuwamura

16 92[28-3-D]

Masato Kimura

16 94[28-4-B]( to be continued )189

Speaker Index
Name

Page no.

Session id

Matthew Min-Hsiung Lin

17 113[28-4-D]

Mayumi Shoji

19 148[29-3-B]

Min-Hung Chen

16 98

[28-4-B]

Mitchell Luskin

11 24[27-1-A]

Mohd Omar

15 75[28-3-C]

Naoki Wada

19 135[29-2-D]

Naoya Yamanaka

15 73[28-3-B]

Norikazu Saito

20 161[29-3-D]

Roberd Saragih

12 41[27-4-B]

Roslinda Mohd Nazar

14 68[28-2-A(2)]Sangho Kum

21 176[29-4-C]

Santiago Madruga

17 106[28-4-C]

Seyoung Park

20 154[29-3-C]

Shin-Hwa Wang

19 146[29-3-B]

Shu-Ming Chang

22 185[29-4-D]

Shugo Manabe

21 172[29-4-B]

Sri Redjeki Pudjaprasetya

17 108[28-4-C]

Swaroop Nandan Bora

13 62

[27-4-D]

Syakila Binti Ahmad

13 56[27-4-D]

Takafumi Miyata

12 34[27-4-B]

Takahito Kashiwabara

20 163[29-3-D]

Takanori Katsura

15 76[28-3-C]

Takashi Ichinomiya

18 126[29-2-B]

Takashi Sakajo

19 141[29-3-B]( to be continued )190

Speaker Index
Name

Page no.

Session id

Takehiko Kinoshita

18 132[29-2-C]

Takeshi Fukaya

12 36[27-4-B]

Takeshi WATANABE

19 150[29-3-B]

Tanay Deshpande

15 72[28-3-B]

Teng-Yao Kuo

16 99[28-4-B]

Tetsuya Ishiwata

20 157[29-3-D]

Tomoyuki Miyaji

19 143[29-3-B]

Tsuyoshi Mizuguchi

16 88[28-3-D]

Victor Didenko

13 58[27-4-D]

Viswanathan Kodakkal Kannan

19 139[29-2-D]

Wah June Leong

20 156[29-3-C]

Wai-Ki CHING

17 117[28-4-D]

Wen-Hann Sheu

11 26[27-2-A(1)]

Wen-Hung Kuo

20 152[29-3-C]

Xiao-Qing Jin

12 38[27-4-B]

Xuefeng LIU

18 133[29-2-C]

Yasumasa Nishiura

14 64[28-1-A]

Yintzer Shih

18 130[29-2-C]

Yongyong Cai

12 28[27-3-A]

Yukihito Suzuki

13 47[27-4-C]

Yung-Ta Li

12 39[27-4-B]

Yusuke Morikura

22 186[29-4-D]

191