Re: [Boost-users] multi-array and pdes

Subject: Re: [Boost-users] multi-array and pdes
From: pmamales_at_[hidden]
Date: 2011-11-07 13:02:06

Hi Larry,
You should really try the paper, not the book (the ppaper is more recent
by 4-5 years).
RK is Runge-Kutta. Basically, one writes the PDE as a system of
ODEs in the open set (w/out the closure). Then advancing over time
is done on this system resulting in the schemes.

My understanding is that the  theta I am referring to is the theta of the numerical stepping, not the
option sensitivity! The statement for second order time accuracy is local! 
According to different thetas you get the global error handling.

What dimension problem are you solving?
What kind of BC's.
It is not clear to me whether one should solve the systems on the open or the closed set.
I am of the opinion that is is the former. It should not matter, if you do Dirichlet BCs.
Make sure though that you employ second order discretization at the boundary (one-sided,of course).

Instead of going with the formal error analysis (always inadequate since a lot of assumptions are tacitly 
imposed) better do the following:
Say you are solving a 2d problem. Write a function u(x,y,t) that you like.
Then derive the PDE and BCs that this satisfy.
Then solve with your machinery.
Compare..
Yes this is nowhere near a mathematical proof but nothing in ADIs really is ;-))
(I know this is worse than others, but will halp you locate the issues: is it BCs? is it oscillations that propagate
because of the one-sided derivatives? you will see).

As far as milti-arrays go, I am very sceptical these days.
This thing looks half=baked..
I do not trust that it can do rotated iterators.
I have resorted to creating the systems to solve on the fly,
by iterating on the proper plane (have created a templated functor that
allows some actor class to do something in the sub-plane that I define).
Typically, the LES(=linear eqtn solver) routines require copying (this is LAPAK for you)
and this seems unavoidable..
 
The tridiagonal systems I am solving either using Thomas algorithm (wiki it), or the LAPACK
tridiagonal solver, which allows for pivoting but is slower. The latter allows for many rhs
vectors, which is useful only if the matrix and the BCs are independent of the dimension(s) on which
you iterate.
My strategy is to solve those over different threads, hence I am stuck with the plumbing.
I am thinking now, if pivoting
I have not looked into your code yet, partially being afraid that it would use C++ features(variadic templates)
that my compiler -unfortunately- does not support.

HTH,
P-
ps: I think this is getting overly specialized for the boost forum. Maybe we should continue privately?
You have my email at home, I think.

---- Larry Evans <cppljevans_at_[hidden]> wrote: 
> On 10/30/11 11:28, petros wrote:
> > Hi Larry,
> > ADI is well suited for problems with sacial and temporal vatiation of
> > the PDE cofficients.
> > Among experts it is indeed considered somewhat of a "hack" (not really
> > though, only not universal solver)
> > and other more elaborate solvers are preferred.
> > Having said all this, it is the industry standard for its speed and
> > ability to tackle a large variety of problems.
> > ADE is a very new method advocated by DD and -to be honest- the fact
> > that noone has picked it up makes me cautius
> > (he-as of now- is the only one supporting it).
> > Would I ever implement it? Sure, after ADI though because of the
> > standard-ness and because once you have one you have
> > enough infra for the other.
> 
> Thanks Petros for that info.
> 
> > Now on the error level you report, it does not make much sense. After
> > all it is a 2nd order scheme, stabilized for time evolution with 2nd
> > order RK.
> 
> What is RK? The only thing I can think of is Runge-Kutta; however, I
> saw nothing about Runge-Kutta in Duffy's chapter 19 on which the ADI
> code I mentioned is based.
> 
> > Are you sure that
> > a) took the BCs into account properly?
> > b) you propagated the BCs as well?
> > c) have you experimented with various values of theta?
> 
> One of the theta's in Duffy's book is one of the Greeks on
> page 131 which is the derivative of option price w.r.t. time.  The
> other might be the theta in equation(7.7) on p. 81; however, that's
> not used in chapter 19.
> 
> Which theta might you be referring to?
> 
> > d) have you experimented with different time step size?
> Yes.
> > HTH
> 
> I've gone over the code again (and again, and again, adding comments
> as I go) trying to figure out what I might have done wrong.  I did
> find that I failed to rotate the delta spaces (corrected now in the
> svn code.  See call to std::rotate); however, I'm still getting 4-6
> times as much error in innocent (the method from section 19.5 of
> Duffy's book) as I'm getting in Duffy's explicit.
> 
> I may just try to understand the books error analysis and try to print
> something about that while the code is running to see why innocent is
> so much worse then explicit; however, that'll take some time :(
> 
> > P-
> > ps: thx for the update. these days I am deep in multithreading LES
> 
> LES = Linear Equation Solver?  The innocent code uses:
> 
> http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/boost/array_stepper/numeric/solve_tridiag.hpp
> 
> which solves a symmetric tridiagonal system.  I did not use mkt or
> ublas or some other library because:
> 
>   1) I wanted to trace the solution and
>      that wouldn't be possible with some library.
> 
>   2) After a brief look at some library(blas or something like that),
>      it looked like it couldn't solve a set of a set of equations, it
>      could only solve 1 set of equations (resulting in a vector
>      result); whereas the solve_tridiag.hpp code can do that,
>      resulting in a multi-dimensional array of results. The number of
>      the sets of equations is the number of times the outermost loop(
>      the one with index, i_node, in triangulate_rhs and
>      back-substitute methods) is executed.  As noted in the comments,
>      I think this outermost loop could be made parallel; hence,
>      multithreading could be used there.
> 
> > (for the purpose of PDE solving) and providing with proper memory
> 
> As I recall, you were planning on using tridiagonal solvers also in
> your PDE solving.  Is that still true?  Are you still planning on
> using multi_array?  If so, are you still needing to rotate the axes?
> If so, how do you do that with multi_array?
> [snip]
> 
> -regards,
> Larry
> 
> 
> 
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users