QuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnQuoteOriginally posted by: outrunIndeed! I see two separate (orthogonal) concepts / libraries* generic numerical derivatives* mesh / grid coordinates mappings / transforms. This will be some mapping between S,t "space" and mesh/grids space i,j (and in doing that also allow for shadow points, 0 or 1 based addressing)The last one has two sub elements: grid, and coordinate transform. For FDM you'll probably have cube-shaped grids, for trees (binomial, trinonimal) you'll have different (triangular) shapes. Both can have the same coordinate transforms. So we have 3 sub libraries:1. numerical derivatives ala dfdx(x1,x2,x3,f1,f2,f3)2. grids, meshes, trees. Basically integer coordinate regions in N dimensions e.g. things like 0<=i<N 3. coordinate tranforms x,y,z <-> f(i,j,k) what do you think?Bang on. All we need now is someone to do itHAHAHAHAHA!!Well, OOglesby sugested 1. We could start a new thread on that, and particullary focus on a good interface. I'm hoping we can provide in general two interfaces: a 'zero-learning-curve' simple one, and a more flexible / generic one. The simple one is basically a synonym for a pre-configured generic out (default policy?)so: *one* pre-cooked interface could be:dfdx(f1,f2,x1,x2) // 2 point, 1st derivativedfdx1(f1,f2,f3,x1,x2,x3) // 3 point, 1st derivative at x1dfdx2(f1,f2,f3,x1,x2,x3) // 3 point, 1st derivative at x2dfdx3(f1,f2,f3,x1,x2,x3) // 3 point, 1st derivative at x3ordfdx3(f1,f2,f3,x1,x2,x3,xi) // 3 point, 1st derivative at xi..don't know.. on a more conceptual level: what are the preconditions? e.g. do we assume x1<x2<x3? Things like that!2. I think we should call these coordinates? Or bound coordinates? 3. I think we should call these coordinate transforms? this one is also a bit tough: I don't think we need bijection? I do think it would be very convenient if we allow for nesting, e.g. use log S (or s/(s+a)) but first apply a cash or stock dividend in S-space. I would also prefer to have the coordinates generic, just specify a good common denominator for an coordinate interface, and then just pick one like double x[3] or std::vector, or boost:tuple but *not* get married to that one (or make divorce possible via traits )I am currently looking at the numerical derivatives task. So far I see several flavors: outrun's pre-cooked interface for 1st and 2nd derivatives an algorithm for N-point, kth derivatives with user specified grid points for x pre-cooked 1st derivative via complex step methods that FinancialAlex has mentioned generalized complex step for N-point, kth derivatives with user specified grid points for xRight now, I am only focusing on flavors 1 and 2.

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

QuoteI'm trying to think how we can fit in Dirichlet, Cauchy, mixed type of boundary conditions. Where would the changes take place, what part in the code?At this stage, my advice would be: don't. The whole issue of boundary conditions is a complicated one and is a separate issue. Both Alan and myself have examined these. It's a long story and has been discussed SSRN and Wilmott July 2011.Are you familiar with the Fichera and Feller approaches? See discussions on Numerical. In some cases, I will show how to 'preprocess' the specific PDE (multiasset, Heston, etc.) so that we deliver what is effectively Dirichlet BC to the bespoke PdeDomain class. A lot of heavy PDE maths.

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: outrunThanks for sharing Daniel! Looks very clean. I'll think a bit about it ok? See if I can shoot at it.* I expect you use the range classes to define upper-lower bounds per axis, e.g t range =[0..2]* Then you'll have something that discretizes it, probably by proving "number of steps" eg N=5 -> dt = 0.5, t[1]=0, t[2]=0.5,...,t[5]=2* And then to pass those values to the boost functionsomerthing = lower_bc( x0, t[ i ]); // for i=1..5Is that how it works more or less (the interaction between range, discretizing, evaluating BC's)?Thijs,You are correct on all fronts. The pde and domain classes are separated for flexiblity. Their instances are created by implemented their boost functions in whatever way you prefer. In my upcoming code I use a namespace containing functions and data along with a factory. On the other side of the pde and domain classes I have various FDM schemes. Very flexible. It can even be used with FEM.I will include provisions for complex coefficients.Now the interesting part. I have many kinds of PDE, so let's take 1 factor problem:u_t = au_xxu_t = a(bu_x)_xu_t = a(bu)_xxu_t = a(b(u)u_x)_x (nonlinear)etc.I do this using Visitor patterns.QuoteFinally, is "Real" a good name that we want to have -- Alan mentioned complex coefficients and multi-precision features of Mathematica, in general we only have different kinds of floating-point-numbers and integers available out of the box. Coefficients can belong to any Ring , including complex, rational and Interval (uncertain coefficients) yes?

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

Yes, agree you need to allow for various non-linearities -- a basic case I would expect QFCL to support is the penalty approach to amer-style options. Other cases of interest are non-linear modelsfor uncertain vol/differential borrowing/lending rates, etc. However, if the non-linearity handling obfuscatesthe basic pde code, maybe that should be an entirely different class of solvers. Ditto for PIDE's.

Last edited by Alan on November 10th, 2011, 11:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: AlanYes, agree you need to allow for various non-linearities -- a basic case I would expect QFCL to support is the penalty approach to amer-style options. Other cases of interest are non-linear modelsfor uncertain vol/differential borrowing/lending rates, etc. However, if the non-linearity handling obfuscatesthe basic pde code, maybe that should be an entirely different class of solvers. Ditto for PIDE's.Clear. These are the modelling issues I am preparing after the initial release. I can model *all* PDEs by using the most general form of coefficients like a(x,t,u, u_x,u_xx) but that will be awful hard to model by FDM. So we need to be realistic and define common categories of linear and NL PDEs, as well as PIDE.A penalty method is a linear PDE + semilinear term: u_t = Lu + f(x,t,u) which I can easily model as a decorator of an existing FDM for linear PDE. Just need to tag on the term f(u). We can reuse code. And throw in BC and domain transformation and then lots of fun starts.Alan,y = S/(S + a) is a good one and this leads to a new PDE, so is y = tanh(aS) leads to a new PDE, ad infinitum.Now, let's say I want a general y = f(S) and S = h(y) and I use these to transform to a *single* PDE(S) to a single PDE(y). Then I would have a PDE that lets user choose own f() and g(). What do you think? It saves a helluva lot of work and manual labour.

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnQuoteI'm trying to think how we can fit in Dirichlet, Cauchy, mixed type of boundary conditions. Where would the changes take place, what part in the code?At this stage, my advice would be: don't. The whole issue of boundary conditions is a complicated one and is a separate issue. Both Alan and myself have examined these. It's a long story and has been discussed SSRN and Wilmott July 2011.Are you familiar with the Fichera and Feller approaches? See discussions on Numerical.No I'm not familiar. I *do* have the Fichera paper though: maybe I should read it.. I ran into these BC when valuing a gas storage (long time ago)Here is Fichera + ADE

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

Another set of scenarios is transforming ('massaging') a PDE in one form into another form, e.g. au_xx + b_ux --> A(Bu_x)_x (integrating factor)A(x)(B(x)u_x)_x --> C(y)(D(y)u_y)_y ( y = x/(x+1))log transform and functional composition of above. Can Heston be put in form u_t = A(Bu_S)_S + C(Du_v)_v + Eu?

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

QuoteOriginally posted by: CuchulainnAlan,y = S/(S + a) is a good one and this leads to a new PDE, so is y = tanh(aS) leads to a new PDE, ad infinitum.Now, let's say I want a general y = f(S) and S = h(y) and I use these to transform to a *single* PDE(S) to a single PDE(y). Then I would have a PDE that lets user choose own f() and g(). What do you think? It saves a helluva lot of work and manual labour.It is tedious to do those sorts of transformations, so helping the user do it would be useful.Reminds me of something Mathematica does, which you might adapt. In Mathematica, if youdo a numerical integration over say (0,xmax) and specificy xmax as some finite number, then it just does itsusual adaptive integration. *But*, if you say xmax=Infinity (Infinity is a Mathematica reserved symbol),then it will automatically do a similar domain transformation on your integrand before handing itoff to the adaptive integrator. [Sometimes this has problems, but it often works fine.]So, for your pde solver, you might consider something similar. Specifying finite boundaries leavesthe problem untouched, but specifying some boundary as +/- Infinity automatically invokes a domaintransform (maybe with an optional Method argument, which specifies which transform, or a default choice).

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

In C++ a templated infinity could be used and it would be in line with the underlying type (e.g. float, double). Polter, is this a good idea?numeric lmits Then we could also model Hull White and log space.

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**60268**Joined:****Location:**Amsterdam-
**Contact:**

QuoteThe Ring is also a good idea.'Ring' could be a concept describing the interface +, *, for mathematical rings (including 0<T> and 1<T> if you get my drift)?? Then real classes implement the concept.

Last edited by Cuchulainn on November 10th, 2011, 11:00 pm, edited 1 time in total.

GZIP: On