A Step-By-Step Procedure To The Numerical Solution For Time-Dependent Partial Derivative Equations in Three Spatial Dimensions
A Step-By-Step Procedure To The Numerical Solution For Time-Dependent Partial Derivative Equations in Three Spatial Dimensions
A Step-By-Step Procedure To The Numerical Solution For Time-Dependent Partial Derivative Equations in Three Spatial Dimensions
Xi Wang
Xinyan Lin
Declaration
We declare that this thesis was composed by ourselves and that the work contained therein is
our own, except where explicitly stated otherwise in the text.
(Xi Wang
Xinyan Lin)
Acknowledgement
To our supervisor Jan R
oman who is the lecturer of Analytical Finance I and II, assistant vice
president and head of market risk in Swedbank in Sweden. First of all, thanks for giving us such
a challengeable topic. Secondly, we appreciate your invaluable experience, patience, encouragement and continuous support. Otherwise, this paper would not have been achieved so far.
Also, thanks to two PHD students: Karl Lundeng
ard and Christopher Engstrom in Malardalen
University, they are also the lecturers in the course Applied Matrix Analysis and Numerical
Method with Matlab. With their experiences and knowledge, the matrices in programming
part have been dealt with successfully.
Abstract
Derivative market has been trading for quite a long time, like options, futures. Meanwhile,
interest rate derivative market is the biggest derivative market among the others in the world.
In recent years, interest rate derivative arises much concern for various purposes. For example,
an investor from one country may have access to getting involved in some derivatives in another
country by exposing himself to the exchange rate. However, in this situation, simulation methods to value such derivative might seem inadequate. On the other hand, numerical methods,
which have done remarkable job when applied to the valuation, therefore, are strongly recommended in the academic community. This paper consolidates insights from another seminar
paper [3] which proposes a Partial Differential Equation (PDE) approach to price cross-currency
interest rate derivatives in terms of three correlated processes, uses finite difference methods
for the spatial discretization of the PDE, and proves that the Alternative Direction Implicit
is the most efficient method to solve such PDE. However, no detailed procedure is shown to
solve the time-dependent parabolic PDE in three spatial dimensions. Thus, our work focuses
on providing step-by-step instructions that enable the reader to follow to solve the parabolic
PDE by finite difference method. Meanwhile, this paper is a guide for readers to see how the
PDE is applicable to general realistic financial problems.
Contents
Abstract
List of Figures
1 Introduction
8
8
8
16
6 Thomas Algorithm
18
7 Summary
19
8 Future Work
20
9 Bibliography
21
A First Appendix
22
B Second Appendix
23
24
D Codes in Matlab
25
List of Figures
3.1
5.1
6.1
Chapter 1
Introduction
In this paper, we investigate the cross-currency swap by using the PDE approach. Particularly,
it is floating to floating cross-currency swap which involves foreign exchange rate and the interest
rates for two currencies. The current standard modeling of such products consist of two onefactor Gaussian models for the term structures and one-factor log-normal model for the spot
FX rate [3][10]. However, these three time-dependent variables construct a three dimension
problem. Generally, such PDE is considered to be a massively computational challenge if it is
solved straightforward. Therefore, we also investigate the efficient and well-known Alternative
Direction Implicit (ADI) method to solve the PDE.
In this paper we avoid to discuss the models and the parameters in the models, but only focus
on the how the PDE approach is applied to find the value of such derivatives. The outline of
the rest of this paper is arranged as follows. In section 2, we present the models and show how
the PDE is derived. Section 3 and 4, Alternative Direction Implicit method and the splitting
method are demonstrated, so that making a bunch of equations into a nice and acceptable
accuracy scheme which results in a convenience computer programming system. Section 5 and
6; generally give a short introduction to the boundary conditions and the Thomas Algorithm.
Section 7 and 8 give a short summary and future work.
Chapter 2
Let u(s, rd , rf , t) denote the domestic value function of a security, where s is spot foreign
exchange rate; rd ,rf are the domestic and foreign short rate respectively [3].
Assume the spot FX as follows [9]:
ds(t) = (rd (t) rf (t))s(t)dt + (t, s(t))s(t)dWs (t)
(2.1.1)
where Ws (t) is a Brownian motion. rd and rf follow the mean-reverting Hull-White model and
have the following dynamics [9]:
drd (t) = (d (t) d (t)rd (t))dt + d (t)dWd (t)
drf (t) = (f (t) f (t)rf (t))dt s ,f (t)f (t)(t, s(t))dt + f (t)dWf (t)
(2.1.2)
(2.1.3)
where deterministic functions i (t) and i (t), with i = d, f , are the mean reversion rates and the
volatility functions. Deterministic function i (t) is a function of time determining the average
direction in which ri moves. dWd (t) and dWf (t) are the other Brownian motions with respect
to rd and rf .
The quanto drift adjustment, s,f f (t)(t, s(t)) for drf (t) comes from changing the measure
from the foreign risk-neutral measure to the domestic risk neutral one. (t, s(t)) is the local
volatility function for the spot FX rate has the functional form [9]
(t, s(t)) = (t)
s(t)
L(t)
(t)1
(2.1.4)
where (t) is the relative volatility function, (t) is the time-dependent constant elasticity of
variance (CEV) parameter and L(t) is a time-dependent scaling constant [3].
2.2
Itos lemma
Because the normalized price process of any security is a martingale under the domestic riskneutral measure, it is easy to apply the Itos lemma and simply set the drift term to be zero.
Itos lemma gives:
du =
+
u
u
u
u
1 2u 2 1 2u 2 1 2u 2
dt +
ds +
ds +
drd +
drf +
dr +
dr
t
s
rd
rf
2 s2
2 rd2 d 2 rf2 f
2u
2u
2u
dsdrd +
dsdrf +
drd drf
rd s
rf s
rf rd
(2.2.1)
u
u
dt +
((rd (t) rf (t))s(t)dt + (t, s(t))s(t)dWs (t))
t
s
u
((d (t) d (t))rd (t)dt + d (t)dWd (t))
rd
u
((f (t) f (t)rf (t))dt f s (t)f (t)(t, s(t))dt + f (t)dWf (t))
rf
1 2u 2
1 2u 2
1 2u 2
2
(
(t,
s(t))s
(t)dt)
+
(
(t)dt)
+
( (t)dt)
d
2
2 s2
2 rd
2 rf2 f
2u
2u
(s,d d (t)(t, s(t))sdt) +
(s,f f (t)(t, s(t))sdt)
srd
srf
2u
(d,f d (t)f (t)dt)
rd rf
(2.2.2)
Rearrange the above equation and group the - dt terms, because of the martingale property, by
setting these terms equal to zero the PDE is derived as:
u
u
u
+ (rd rf )s
+ (d (t) d (t)rd )
t
s
rd
u
(f (t) f (t)rf f,s f (t)(t, s(t)))
rf
1 2
2u 1
2u
2u 1
(t, s(t))s2 2 + d2 (t) 2 + f2 (t) 2
2
s
2
rd
2
rf
rd u
2u
2u
2u
+ s,f f (t)(t, s(t))s
+ d,f d (t)f (t)
srd
srf
rd rf
(2.2.3)
Chapter 3
finite differences
In the discretization, the first and second derivative terms can all be approximated by the
central difference method as presented following [3][7]:
Similarly,
m
m
m
2
um
um
u
i+1,j,k ui1,j,k u
i+1,j,k 2ui,j,k + ui1,j,k
, 2
s
24s
s
4s2
(3.1.1)
m
m
m
2
um
um
u
i,j+1,k ui,j1,k u
i,j+1,k 2ui,j,k + ui,j1,k
, 2
rd
24rd
rd
4rd2
(3.1.2)
m
m
m
2
um
um
u
i,j,k+1 2ui,j,k + ui,j,k1
i,j,k+1 ui,j,k1 u
, 2
rf
24rf
rf
4rf2
(3.1.3)
m1
um
u
i,j,k ui,j,k
t
4t
(3.1.4)
(3.1.5)
m
m
m
um
2u
i+1,j,k+1 + ui1,j,k1 ui1,j,k+1 ui+1,j,k1
=
rf s
44s4rf
(3.1.6)
m
m
m
um
2u
i,j+1,k+1 + ui,j1,k1 ui,j1,k+1 ui,j+1,k1
=
rd rf
44rf 4rd
(3.1.7)
10
um
i,j1,k1
i,j1,k+1 ui,j+1,k1
+ d,f d (t)f (t)
44rf 4rd
rd u
!
m1
um
i,j,k ui,j,k
=
(3.1.8)
4
Rearrange in terms of the directions i , j , k :
(rd rf )s 2 (t, s(t))s2
(rd rf )s 2 (t, s(t))s2
m
+
+
u
+
um
i+1,j,k
i1,j,k
24s
24s2
24s
24s2
d (t) d (t))rd
2 (t)
2 (t, s(t))s2 m
ui,j,k +
+ d 2 um
i,j+1,k
2
4s
24rd
24rd
d (t) d (t))rd
d2 (t) m
2 (t)
+
u
+ d 2 um
i,j1,k
24rd
24rd
4rd2 i,j,k
!
f2 (t) m
f2 (t)
f (t) f (t)rf f,s f (t)(t, s(t))
m
+
u
u
+
i,j,k+1
24rf
24rf2
4rf2 i,j,k
!
f2 (t)
f (t) f (t)rf f,s f (t)(t, s(t))
um
+
+
i,j,k1
24rf
24rf2
s,d d (t)(t, s(t))s m
m
m
ui+1,j+1,k + um
i1,j1,k ui1,j+1,k ui+1,j1,k
44s4rd
s,f f (t)(t, s(t))s m
m
m
ui+1,j,k+1 + um
+
i1,j,k1 ui1,j,k+1 ui+1,j,k1
44s4rf
d,f d (t)f (t) m
m
m
+
ui,j+1,k+1 + um
i,j1,k1 ui,j1,k+1 ui,j+1,k1
44rf 4rd
rd u
m1
um
i,j,k ui,j,k
=
4
+
3.2
(3.1.9)
The basic idea behind the ADI method is to split the finite difference equations into different
directions. Since in our case, the three dimensions form a cube for every time step, Figure 3.1
demonstrates the idea.
For each split, only one direction is taken implicit. As the procedure continues, the exact
solution will be found. From the finite difference method, Equation 3.1.9 can be written in the
11
um um1
4
(3.2.1)
m
m
m
Where Am
1 ,A2 ,A3 denote the terms in each direction. A0 denotes the cross derivative terms,
m
m
are treated in an explicit fashion. And assume the term rd u is distributed evenly over Am
1 ,A2 ,A3
[3]. Therefore:
u 2 u 1
+ 2 rd
s
s
3
u
2u 1
+ 2 rd
rd
rd
3
Am
1
Am
2
Am
3
u
2u 1
+ 2 rd
rf
rf
3
Am
0
2u
2u
2u
+
+
rd s rd rf
rf s
(3.2.2)
2 s2
um
i,j,k
4s2
d2
m
= c um
um
i,j,k
i,j+1,k + d ui,j1,k
4rd2
m
= a um
i+1,j,k + b ui1,j,k
m
= e um
i,j,k+1 + f ui,j,k1
1
rd um
i,j,k
3
1
rd um
i,j,k
3
f2
1
m
um
i,j,k rd ui,j,k
4rf2
3
m
m
In such form, the following matrices represent the Am
1 ,A2 ,A3 respectively:
2 2
4ss2 13 rd a
0
..
..
Am
.
1 =
b
.
2 2
0
4ss2 13 rd
Am
2 =
4rd2 13 rd
d
0
..
0
..
.
d2
4r2 13 rd
d
12
(3.2.3)
4rf2 31 rd
..
Am
3 =
0
..
.
.
2
4rf2 13 rd
f
Where
a
b =
c =
(rd rf )s
2 s2
+
24s
24s2
(rd rf )s
2 s2
+
24s
24s2
(d d rd )
d2
+
24rd
24rd2
(d d rd )
d2
+
24rd
24rd2
d =
e =
f2
(f f rf s,f f )
+
24rf
24rf2
f2
(f f rf s,f f )
+
24rf
24rf2
13
Chapter 4
The HV scheme
The Hundsdorfer and Verwer (HV) [5] splitting scheme as suggested in [3] approximates the
exact solution U m :
V0
4
Am
(I
i )Vi
2
V0
(I 12 4 Am
i )Vi
m
U
=
=
=
=
=
U m1 + 4(Am1 U m1 + g m1 )
Vi1 21 4 Am1
+ 12 4 (gim gim1 ),
i = 1, 2, 3
i
1
m
V0 + 2 4 (A V3 Am1 U m1 ) + 12 4 (g m g m1 )
Vi1 21 4 Am
i = 1, 2, 3
i V3 ,
V3 .
Here the vector g m is the boundary conditions for the corresponding directions, will be discussed
in the later section. In our case, the boundary conditions form a cube. Besides, from now on,
the A matrices are defined in a slightly different manner because of the computational and
programming purpose; they are presented in the Appendix A.
4.1
Applying the Douglas [4] ADI scheme to the equation 3.2.1, obtain a first evaluate of the
solution at time m + 1 by taking implicitly half of the A1 :
1
1
um+ 3 um
1
A1 (um+ 3 + um ) + A2 um + A3 um + A0 um =
2
4
(4.1.1)
1
2
1
1
um+ 3 um
A1 (um+ 3 + um ) + A2 (um+ 3 + um ) + A3 um + A0 um =
2
2
4
(4.1.2)
1
2
1
1
1
um+1 um
A1 (um+ 3 + um ) + A2 (um+ 3 + um ) + A3 (um+1 + um ) + A0 um =
2
2
2
4
(4.1.3)
The intermediate values um+ 3 , um+ 3 can be eliminated. The reader is strongly recommended
to the reference [4] for detailed elimination procedures.
For splitting purpose, the above can be rewritten as:
1
1
(I 4 A1 )um+ 3
2
2
1
(I 4 A2 )um+ 3
2
1
(I 4 A3 )um+1
2
=
=
=
1
um + 4 (A2 um + A3 um + A0 um ) + 4 A1 um
(4.1.4)
2
1
1
1
um + 4 A1 (um+ 3 + um ) + A2 um + A3 um + A0 um (4.1.5)
2
2
1
2
1
1
um + 4 A1 (um+ 3 + um ) + A2 (um+ 3 + um )
2
2
1
+ A3 um + A0 um
(4.1.6)
2
14
(4.1.7)
1
m+ 2
m+ 2
m+1 m+1
3
3
ui,j1,k
,ui,j,k3 ,ui,j+1,k
,um+1
i,j,k1 ,ui,j,k ,ui,j,k+1 are unknown. As stated earlier, A1 , A2 , A3 are
tri-diagonal matrices, therefore, the equations 4.1.4, 4.1.5, 4.1.6 can be expressed below:
m+ 1
u1,j,k3
m+ 13
u
2,j,k
m+ 13
u3,j,k
1
I 4 A1
um+ 13
2
4,j,k
..
m+ 13
ui,j,k
m+ 2
ui,1,k3
m+ 23
u
i,2,k
m+ 23
ui,3,k
1
I 4 A2
um+ 23
2
i,4,k
..
m+ 23
ui,j,k
m+1
ui,j,1
um+1
i,j,2
um+1
i,j,3
1
I 4 A3
um+1
2
i,j,4
..
.
um+1
i,j,k
= [1 ]
= [2 ]
= [3 ]
15
Chapter 5
= g1
s
24s
m
m
m
ui2,j,k 4 ui1,j,k + 3 ui,j,k
u
= g2
s
24s
f orward
backward
m
m
m
The start and end surfaces are um
0,j,k , uI,j,k and then, the u0,j,k and uI,j,k are solved from the
forward and backward differences which will result in the different coefficients in the matrix on
the left hand side of the equation as well as the right hand side. This boundary condition is
more complicated in implementations. We will continue to discuss this topic in the next paper,
where we will do the calibrations and show different boundary conditions result in comparable
results. For now, in this paper, we simply apply the Dirichlet boundary conditions where there
are fixed values on the end surfaces. This type of the boundary condition takes the values on
the boundary points as known, but only changes the right hand side of the equation when the
boundary conditions are considered.
F1 g1
F2
F
3
..
..
Fn1
Fn gn
To be more explicit, in our case, the boundary conditions are surfaces. Figure 5.1 below
illustrates the boundaries and the area we are really concern about which is indicated by the
red hypercube. Readers can get helps to understand the boundary conditions by having a look
at the [1, 2, 6] in the reference.
16
Figure 5.1: Illustration on the boundary conditions and the solution area
17
Chapter 6
Thomas Algorithm
HV scheme has made the algorithm clear enough to solve the U m , but since the algorithm
involves solving a large number of systems of equations which may cause difficulties in the
implementation. Fortunately, knowing that A1 , A2 , A3 are three tri-diagonal matrices, the
Thomas algorithm [11] has been invented and recommended to help us to do the job. The basic
idea behind the algorithm is substitution which will result in a new tri-diagonal matrix. Of
course, the vector on the right hand side of the equation will be substituted as well. However,
the advantage is, the elements in sub-diagonal in new tri-diagonal matrix are zeros, and the
diagonal elements are ones. Figure 6.1 illustrates how the diagonal matrix is transformed using
the Thomas algorithm.
18
Chapter 7
Summary
Up until here we have made complete procedures to the numerical method for three dimensional
partial derivative equations. As can be seen, such problems could be extremely complex. The
matrices we illustrated on this paper only reflect the general forms. When the indices are
considered, we are actually dealing with hundreds of equations. Thus, to find the final exact
solution, every step is critical; our goal is to efficiently obtain the solution with high accuracy.
Thanks to the other academicians; the problems in the real situations can successfully be solved
even with different order of errors. The programming languages have been playing a vital part,
for example the MATLAB, C++, etc., in the appendix; we enclose the MATLAB program for
solving such problems. However, the program has to be modified for any other specific type of
instruments. Therefore, in this paper, we skip analyzing the results from the program but only
give a short insight into it. More details will be covered in the next paper.
19
Chapter 8
Future Work
So far, this paper only gives the general understanding to the numerical method of solving partial
differential equations in three dimensions. As has been said in the introduction, this paper has
not talked about any parameters in the formulas; obviously, those parameters will significantly
influence the results. In the next paper, will test the program and try to calibrate the parameters
from the real market data. When dealing with the real market data, the boundary conditions
we consider in the paper may not be suitable anymore because of the time-dependent variables.
Thus, to achieve great accuracy, more specific boundary conditions will be included. Finally,
the MATLAB program built in the way that suits for general instruments whichever involves
three dimensions, estimating some of such popularly traded instruments will be illustrated in
the next paper as well.
20
Chapter 9
Bibliography
[1] Nasser M. Abbasi. Neumann boundary conditions on 2d grid with non-uniform mesh space
for elliptic pde, 2012.
[2] Wolfgang Arendt and Mahamadi Warma. Dirichlet and neumann boundary conditions:what is in between? Journal of Evolution Equations, 3:119135, 2003.
[3] Duy Minh Dang, Christina C. Christara, Kenneth R. Jackson, and Asif Lakhany. A pde
pricing framework for cross-currency interest rate derivatives. Procedia CS, 1(1):23712380,
2010.
[4] Jim Douglas. Alternating direction methods for three space variables. Numerische Mathematik, 4:4163, 1962.
[5] K.J. in t Hout and B.D. Welfert. Unconditional stability of second-order adi schemes
applied to multi-dimensional diffusion equations with mixed derivative terms. Applied
Numerical Mathematics, 59(34):677 692, 2009. ce:titleSelected Papers from NUMDIFF11/ce:title.
[6] J.Izadian and S.S.Jalalian. A new method for solving 3d elliptic problem with dirichlet
or neumann boundary conditions using finite difference method. Applied Mathematical
Sciences, 6(34):16551666, 2012.
[7] Zhilin Li. Finite difference method basics.
[8] Sensen Lin. Finite difference schemes for heston model, 2008.
[9] Vladimir Piterbarg. Smiling hybrids. Risk magazine, pages 6670, May 2006.
[10] Jason Sippel and Shoichi Ohkoshi. All power to prdc notes. Risk magazine, pages 3133,
November 2002.
[11] W.T.Lee. Tridiagonal matrices: Thomas algorithm.
21
Appendix A
First Appendix
For the programming purpose, the following transformations are needed which are based on
Lins work[8]. We let:
4t
4s
4t
m
m
m
= Rs , 4s
2 = Rss , ui+1,j,k ui1,j,k = s ui,j,k ,
m
m
m
um
i+1,j,k 2ui,j,k + ui1,j,k = ss ui,j,k
4t
4rd
4t
m
m
m
= Rd , 4r
2 = Rdd , ui,j+1,k ui,j1,k = d ui,j,k ,
d
m
m
m
um
i,j+1,k 2ui,j,k + ui,j1,k = dd ui,j,k
4t
4rf
4t
m
m
m
= Rf , 4r
2 = Rf f , ui,j,k+1 ui,j,k1 = f ui,j,k ,
f
m
m
m
um
i,j,k+1 2ui,j,k + ui,j,k1 = f f ui,j,k
4t
4t
4t
4rd 4s = Rsd , 4rf 4s = Rsf , 4rd 4f = Rdf
m
m
m
m
um
i+1,j+1,k + ui1,j1,k ui1,j+1,k ui+1,j1,k = sd ui,j,k
m
m
m
m
um
i+1,j,k+1 + ui1,j,k1 ui1,j,k+1 ui+1,j,k1 = sf ui,j,k
m
m
m
m
um
i,j+1,k+1 + ui,j1,k1 ui,j+1,k1 ui,j1,k+1 = df ui,j,k
A1 =
A2 =
2 s2
1
2 Rss ss 3 rd
2
d d rd
Rd d + 2d Rdd dd 13 rd
2
(rd rf )s
Rs s
2
2
f f rf s,f f
Rf f + 2f Rf f f f 31 rd
2
s,d d s
s
Rs ds d + s,f 4 f Rsf sf + d,f 4 f d Rdf df
4
A3 =
A0 =
2 2
(rd rf )s
2 s2
4ss2 13 rd
+ 24s
0
2
24s
..
..
rf )s
2 s2
.
A1 = (rd24s
+ 24s
.
..
2 2
.
4ss2 13 rd
d2
2
d d rd
4rd2 13 rd
+
0
24rd
24rd2
d
r
2
..
.
.
d
d d
d
.
A2 =
+
.
24rd
24rd2
..
d2
.
4r2 31 rd
d
f2
f2
f f rf s,f f
1
4r2 3 rd
+ 24r
0
2
24r
f
f
f
f f rf s,f f
..
..
f2
.
A3 =
+ 24r
.
2
24rf
f
2
..
4rf2 13 rd
f
22
Appendix B
Second Appendix
For the boundary conditions, approximations to the partial derivative have to be forward and
backward second-order one-sided finite difference. More specifically, for the boundary on the
first points, use the forward difference, for the last points, use the backward points Secondorder
onesided approximation for the derivatives:
m
m
4 um
u
i+1,j,k ui+2,j,k 3 ui,j,k
s
24s
m
m
um
u
u
i2,j,k
i1,j,k + 3 ui,j,k
s
24s
m
m
4 um
u
u
i,j+1,k
i,j+2,k 3 ui,j,k
rd
24rd
m
m
um
u
i,j2,k 4 ui,j1,k + 3 ui,j,k
rd
24rd
m
m
4 ui,j,k+1 um
u
i,j,k+2 3 ui,j,k
rf
24rf
m
m
ui,j,k2 4 um
u
i,j,k1 + 3 ui,j,k
rf
24rf
23
f orward
backward
f orward
backward
f orward
backward
Appendix C
24
Appendix D
Codes in MATLAB
function f = TransformA0 ( u0 , rhosd , r h o s f , r h o d f , sigmad , s i g m a f ,gamma, dt )
clear f
ns = s i z e ( u0 , 1 ) 1 ;
nrd = s i z e ( u0 , 2 ) 1 ;
n r f = s i z e ( u0 , 3 ) 1 ;
s l i n e = 0 : 1 : ns ;
d l i n e = 0 : 1 : nrd ;
f l i n e = 0:1: nrf ;
u1 = u0 ( 3 : ns + 1 , 3 : nrd + 1 , 2 : n r f ) . . .
+ u0 ( 1 : ns 1 , 1 : nrd 1 , 2 : n r f ) . . .
u0 ( 1 : ns 1 , 3 : nrd + 1 , 2 : n r f ) . . .
u0 ( 3 : ns + 1 , 1 : nrd 1 , 2 : n r f ) ;
u2 = u0 ( 3 : ns + 1 , 2 : nrd , 3 : n r f + 1 ) . . .
+ u0 ( 1 : ns 1 , 2 : nrd , 1 : n r f 1 ) . . .
u0 ( 1 : ns 1 , 2 : nrd , 3 : n r f + 1 ) . . .
u0 ( 3 : ns + 1 , 2 : nrd , 1 : n r f 1 ) ;
u3 = u0 ( 2 : ns , 3 : nrd + 1 , 3 : n r f + 1 ) . . .
+ u0 ( 2 : ns , 1 : nrd 1 , 1 : n r f 1 ) . . .
u0 ( 2 : ns , 1 : nrd 1 , 3 : n r f + 1 ) . . .
u0 ( 2 : ns , 3 : nrd + 1 , 1 : n r f 1 ) ;
u4 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
u4 ( 2 : ns , 2 : nrd , 2 : n r f ) = u1 ;
u5 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
u5 ( 2 : ns , 2 : nrd , 2 : n r f ) = u2 ;
u6 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
u6 ( 2 : ns , 2 : nrd , 2 : n r f ) = u3 ;
sd = s l i n e d l i n e ;
sf = sline f l i n e ;
df = dline f l i n e ;
f 1 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
f 2 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
f 3 = zeros ( ns + 1 , nrd + 1 , n r f + 1 ) ;
f o r i = 1 : ns + 1
f 1 ( : , : , i ) = sd u4 ( : , : , i ) ;
f 2 ( : , : , i ) = s f u5 ( : , : , i ) ;
f 3 ( : , : , i ) = d f u6 ( : , : , i ) ;
end
f 1 1 = r h o s d sigmad gamma dt f 1 / 4 ;
f 2 2 = r h o s f s i g m a f gamma dt f 2 / 4 ;
f 3 3 = r h o d f s i g m a f sigmad dt f 3 / 4 ;
f0 = f11 + f22 + f33 ;
f = f0 ;
end
25
s i z e ( u0 , 1 ) 1 ;
s i z e ( u0 , 2 ) 1 ;
s i z e ( u0 , 3 ) 1 ;
= 1 : ns 1;
f = zeros ( s i z e ( u0 ) ) ;
B1 = s l i n e . 2 gamma. 2 / 2 ;
B1 = [ B1 ; 0 ; 0 ] ;
B2 = s l i n e . 2 gamma. 2 ;
B2 = [ 0 ; B2 ; 0 ] ;
B3 = s l i n e . 2 gamma. 2 / 2 ;
B3 = [ 0 ; 0 ; B3 ] ;
B = spdiags ( [ B1 B2 B3 ] , [ 1 0 1 ] , ns + 1 , ns + 1 ) ;
C1
C1
C2
C2
C3
C3
C
=
=
=
=
=
=
=
( rd r f ) s l i n e / 2 ;
[ C1 ; 0 ; 0 ] ;
rd o n e s ( ns 1 , 1 ) / 3 ;
[ 0 ; C2 ; 0 ] ;
( rd r f ) s l i n e / 2 ;
[ 0 ; 0 ; C3 ] ;
spdiags ( [ C1 C2 C3 ] , [ 1 0 1 ] , ns + 1 , ns + 1 ) ;
for j = 2
for k =
P1
G1
temp1
f (2 :
end
end
end
: nrd
2 : nrf
= u0 ( : , j , k ) ;
= B + C;
= dt G1P1 ;
ns , j , k ) = temp1 ( 2 : ns ) ;
=
=
=
=
=
=
=
for i = 2
for k =
P2
temp2
:
2
=
=
ns
: nrf
u0 ( i , : , k ) ;
dt DP2 ;
26
: ns
2 : nrd
= u0 ( i , j , : ) ;
= reshape ( P3 , [ n r f + 1 , 1 ] ) ;
= dt EP3 ;
j , 2 : n r f ) = temp3 ( 2 : n r f ) ;
27
28
29