Serving the Quantitative Finance Community

Arnheim
Posts: 42
Joined: July 11th, 2003, 10:36 am

Finite difference - CFD technique

Helloa) Dirichlet Boundary is straightforward. More tricky is linearity.Managed to get Quasi-Linearity running (in the sense of Tavella/Randall), ie discretize PDE at boundary, skip convex terms.(Explicit before you do sweep, Implicit in resp dimension after the sweep for the opposite boundaries (we cant go to the boundary itself with ADE)).This works reasonably well.Does anybody have an idea how to do u_xx=0 eg?b) ADE 2D with explicit cross term is not uncond stable anymore.Got some situations where it blew, small space step compared to time step. Generally no problem.c) I find it necessary to do Extrapolation. Scheme only does not yield desired accuracy. Is there a way to improve on spatial convergence as well (maybe similar to Richardson in time)?d) Found below discretization for cross term with mixed impl/expl terms.Is anybody using this for ADE? Cant see how, we have i+1 or j+1, so we cant be expl anymoreWill try 3D nowRgds

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

QuoteOriginally posted by: ArnheimHelloa) Dirichlet Boundary is straightforward. More tricky is linearity.Managed to get Quasi-Linearity running (in the sense of Tavella/Randall), ie discretize PDE at boundary, skip convex terms.(Explicit before you do sweep, Implicit in resp dimension after the sweep for the opposite boundaries (we cant go to the boundary itself with ADE)).This works reasonably well.Does anybody have an idea how to do u_xx=0 eg?b) ADE 2D with explicit cross term is not uncond stable anymore.Got some situations where it blew, small space step compared to time step. Generally no problem.c) I find it necessary to do Extrapolation. Scheme only does not yield desired accuracy. Is there a way to improve on spatial convergence as well (maybe similar to Richardson in time)?d) Found below discretization for cross term with mixed impl/expl terms.Is anybody using this for ADE? Cant see how, we have i+1 or j+1, so we cant be expl anymoreWill try 3D nowRgdsThanks, congratulations you are now the ADE expert (I a just a simple trainer).Some answers based on stuff from Sammus and DJD:In the meantime I have worked out the UXX = 0 at boundary by a kind of perturbation because we have the terms U(n+1, 1) and U(n+1, 0) which destroys explicit. HOWEVER, I perturb U(n+1, 1) to something in U(n,1) so only U(n+1, 0) is left and solvable on LHS, so back to old Dirichlet again. The resulting matrix is also an M-matrix.Also done Ut + aUx using the same ANE trick and have proved stability using von Neumann and Positivity stuff as you know. Next step is to do it for Ut + aUx + bUy but in a few weeks (in London next week). Sammus needs this for Heston.I have a design pattern framework and will plug in the ANE stuff soon. Further, I want to do a theoretical numerical analysis as well of the problem.For the cross term I think we have the same conclusions as yourself (a combination of IMEX)?Does this mean the Death of ADI? Sammus and myself are possibly only ones working on ANE (ADE). This is obscure work from the Siberian hinterland of the 60's and see Saul'yevs book. It is used in real CFD.We are interested in exchange of results. Viele gruessen.P.S. Have not got to extrapolation but I think that it is possible after we get the O(h) stuff.
Last edited by Cuchulainn on March 8th, 2005, 11:00 pm, edited 1 time in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

sammus
Posts: 319
Joined: November 11th, 2003, 6:21 am

Finite difference - CFD technique

Hi yall ,Just got back from work. The original ADE is introduced by Saul'yev as such: , where the term is replaced by . Therefore, I guess the first order could be written as , then we switch around the explicit and implicit terms, or maybe ?? In a similar manner, the cross term could be thought of as , where we have both explicit and implicit terms. Hopefully such scheme is unconditionally stable. We will see. For the 3D prob, I doubt the ADE method is applicable since we cannot use the central difference there. yes no yes no....
Last edited by sammus on March 9th, 2005, 11:00 pm, edited 1 time in total.

sammus
Posts: 319
Joined: November 11th, 2003, 6:21 am

Finite difference - CFD technique

I had some thoughts while I was trying to sell alcohol to an under-age. I must've been crazy. alright, for a parabolic PDE of the Heston type , the ADE method indicates that for the first sweep we have ^ * *^ * *^ ^ * which requires the computation from the top to the bottom and from the right to the left. We start with solving the point at the upper-right corner (we define i goes from 1 to N and j goes from 1 to M) and then we go down before we move on to the left column. If we introduce the fictitious points, then the middle point can be determined by 8 surrounding points plus its own value at the previous time level. Ignore the 4 explicit points since they are not a problem. If we discretize the zero convexity condition and the same way those implicit and explicit terms are defined in the main equation, the ghost points and can be obsorbed by known interior points and itself. Then the "corner"ghost point U_(N+1,M+1)^(n+1) can also be eliminated. In order to kill off the last ghost point U_(N+1,M-1)^(n+1), we need to know the value of the point U_(N,M-1)^(n+1) which is to be solved next step. Therefore, we build a bi-diagonal matrix to solve for the rightmost boundary points implicitly. After getting the boundary values, we move to the left and the explicit life is always good. For the second sweep, we have* ^ ^* * ^* * ^Similarly, we start from the bottom to the top and from the right to the left. For the Heston model, we have Dirichlet at the bottom, ie. V=0 when S=0. We take natural conditions for vol=0. Fortunetely, the PDE for that condition doesn't come with second and mixed derivatives, so we have ()^ * * ^()*at the leftmost boundary(() means nothing there). Take the starting point U_(1,2)^(n+1) as an example. We can let the "corner" ghost point U_(0,1)^(n+1) be zero in consistency with the Dirichlet condition. Another ghost point U_(0,2)^(n+1) can be completely absorbed according to the natural condition. For the last ghost point U_(0,3)^(n+1), the above graph tells us that we need other 4 interior points. We can again solve the leftmost boundary values from the bottom to the top by constructing another bi-diagonal matrix. After that, the rest of the scheme would be totally explicit. Those are just my rough ideas. We probably need to define the initial ghost point values as a part of initial conditions. With directional computation we can avoid solving some crazy matrices. I think the ADE for the Heston model is workable but I have no clue of its robustness. Hope I made myself clear. Any comments and hints are very much welcome. The Latex editor is not working at this moment. Sorry for the sloppy notations. I will fix them later. gotta go sleep. see you guys.
Last edited by sammus on March 10th, 2005, 11:00 pm, edited 1 time in total.

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

Sammus, Arnheim,good dayHere are some thoughts, answers, dos and donts, questions preliminary results on ADE/ANE. I have done a bit of numerical analysis in this stuff in the last 2 weeks, but not enough maybe to be coherent. My snippets might be for further discussion. Let’s focus on 1 and 2 factors for moment, OK? And let’s for moment forget about cross derivatives.My Facts (are they true)1. ADE works OK for Dirichlet BC because we are sweeping in the Saul’yev method.2. ADE is a bit low order O(h + k), but it works. Later we can do extrapolation.3. The first term Ut + aUx + bUy is a combination of IMEX.4. Convexity BC Uxx = 0 via ghost points is OK albeit that an ‘interior’ implicit value U(n+1, 1) is introducedWhat have I done5. Offload process: the value U(n+1, 1) is perturbed to U(n,1) so that we have to solve for U(n,0) on boundary. I have done a stability analysis and looks OK6. Have not programmed it (yet) in C++7. I have applied Saul’yev to Ut +aUx (!!!!! This is a wave equation, FDM is tricky) and stability is OK. Must do it for Heston in 2 factorsProposalWe should stick to pure ADE for moment. If it does not work we might try a weak IM with IM on the boundary and EX in interior as Sammus suggests.Unfor. Will not be able to work on this next week because of training. Will look at posts however.P.S. Anyone have experience with these FDM issues?P.P.S. using centred differnces j-1, j+1 for first order PDE on biundary is tricky IMHO. I would use one-sided O(h) for moment. Be aware of characteristic direction like what Arnheim and I discussed a while ago.
Last edited by Cuchulainn on March 11th, 2005, 11:00 pm, edited 1 time in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

Arnheim, SammusI have done a Saul'yev on the initial Boun val prob Ut + aUx = 0, a> 0U(0, t) = g(t) // BC FROM LEFT HAND SIDE!!!!!using bc from leftInstead of FTCS, i.e.(U(n+1, j) - U(n,j))/k + a (U(n+1, j+1) - U(n+1, j-1))/(2h)I use (U(n+1, j) - U(n,j))/k + a (U(n, j+1) - U(n+1, j-1))/(2h), j >=1 and have proved that its VNeumann stab. fac. is 1. I think you might need this for a) first order term in BSb) BC for Heston It is a good idea I think to look at this kind of problem on its own merits because this is where we are having problems. keep plugging.
Last edited by Cuchulainn on March 11th, 2005, 11:00 pm, edited 1 time in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

Sammus,I have read Roache "Fund. of CFD" 1998 and he says:1. ADE for Dirichlet is explicit but for other BC we have implicit at boiundary2. Nobody has done ADE well for convetive-diffusionAt some cut-off date we might have to think about 1 if the "fully" explicit does not work. There are search libraries for finding articles related to above?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

sammus
Posts: 319
Joined: November 11th, 2003, 6:21 am

Finite difference - CFD technique

Hi Daniel, Thanks for your excellent explanation. I will do some more readings about ADE and implement the above schemes to see what the results will be. Hoepfully we can do something on the convective-diffusion in ANE/ADE. That is stimulating. Cheers,
Last edited by sammus on March 14th, 2005, 11:00 pm, edited 1 time in total.

sammus
Posts: 319
Joined: November 11th, 2003, 6:21 am

Finite difference - CFD technique

>I use (U(n+1, j) - U(n,j))/k + a (U(n, j+1) - U(n+1, j-1))/(2h), j >=1 also looks like one sweep of ADE at the boundary to me. Should we do one more sweep the other way around? It may be trivial.

Arnheim
Posts: 42
Joined: July 11th, 2003, 10:36 am

Finite difference - CFD technique

Ad ADE:1. Linearity at boundaryI tried linearity (u_xx=0) in the point just before boundary. Then we have two conditions on point x_1: One from the PDE in u_1 and u_0 using ADE Operators for second order deriv. The other from the linearity condition in x_1! This gives us two equations in u_0 and u_1 and can be solved expl.This idea also works in multidimension. We have to go along the boundary first and solve those equations explicitly. Little system in the corner (dimension+1).Works fairly well.Didnt quite understand how you formulate linearity condition. Could you pin down the math?2. Reaction termTried theta-scheme (mixed impl-expl in time) on reaction term. Implicit reaction term consitently gives best result. Scheme still expl.3. Advection termThats a tricky one. Tried the following variations on the first order term:a) Explicit one-sided in spaceb) Explicit central in spacec) ADE kind operator as you wrote it downd) Centered in time, centered in spacee) centered in time, forward in spacef) centered in time, backward in spaceIt turns out, that expl works best (a&b) for a wide range of advection/diffusion equations. And in fact, c) is worst. Behavior very dependent whether PDE is advection or diffusion dominated. Dont have a definite conclusion yet.Would be interested if you have similar results. Will work explicit in advection for the time being.Rgds

Arnheim
Posts: 42
Joined: July 11th, 2003, 10:36 am

Finite difference - CFD technique

Another observation ad advection:If d/dx>>0 explicit seems to work onlyIf d/dx<<0 implicit or ADE operator work betterDoes anybody have an explanation?

Arnheim
Posts: 42
Joined: July 11th, 2003, 10:36 am

Finite difference - CFD technique

Another observation ad advection:If d/dx>>0 explicit seems to work onlyIf d/dx<<0 implicit or ADE operator work betterDoes anybody have an explanation?

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

QuoteOriginally posted by: ArnheimAnother observation ad advection:If d/dx>>0 explicit seems to work onlyIf d/dx<<0 implicit or ADE operator work betterDoes anybody have an explanation?HII do not understand the notation d/dx. You mean the coefficient of the first derivative? Has it to do with upwinding stuff? Does << mean <=?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

3. Advection termThats a tricky one. Tried the following variations on the first order term:a) Explicit one-sided in spaceb) Explicit central in spacec) ADE kind operator as you wrote it downd) Centered in time, centered in spacee) centered in time, forward in spacef) centered in time, backward in spaceIt turns out, that expl works best (a&b) for a wide range of advection/diffusion equations. And in fact, c) is worst. Behavior very dependent whether PDE is advection or diffusion dominated. Dont have a definite conclusion yet.Would be interested if you have similar results. Will work explicit in advection for the time being.Ok, one less option to look at. Is my scheme stable? I did a von Neumann and was neutrally stable (did not program it yet, however).Maybe explicit for moment.I don't propose it just yet but exponentail fittting FDM is good for advection or diffusion dominated. Might be a sledgehammer, solution should be easier.P.S> I would be inclined to do one-step in both x and t for monent until the path is clear.------------------
Last edited by Cuchulainn on March 13th, 2005, 11:00 pm, edited 1 time in total.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl

Cuchulainn
Posts: 64980
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Finite difference - CFD technique

Sammusdifferent questionImportantSomeone want to program ADI for 2 factor. 3 good reasons for doing/not doing it.Everyone says ADI, hello I do ADI, have you seen the article on ADI, I am going to do an ADI ...In my experience it is inferior to other FDMAny punch lines?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl