Riccati ODEs

Ramiro Vignolo

unread,

Apr 23, 2019, 3:36:16 PM4/23/19

to was...@seamplex.com

Hi all! It's been a while...

As you may recall, sometimes I come around and ask for your help in some kind of, what I think, is a nice problem. The last time was about a puzzle (which I have solved XD).

Let's focus on today's problem. I have to solve the following system of Riccati ODEs:

$\frac{{{d}{A}{\left({t}\right)}}}{{{\left.{d}{t}\right.}}}-\kappa{\left({t}\right)}\theta{\left({t}\right)}{B}{\left({t}\right)}+\frac{1}{2}{\sigma}^{2}{\left({t}\right)}\alpha{\left({t}\right)}{B}{\left({t}\right)}={0},$

$-\frac{{{d}{B}{\left({t}\right)}}}{{{\left.{d}{t}\right.}}}-\kappa{\left({t}\right)}{B}{\left({t}\right)}+\frac{1}{2}{\sigma}^{2}{\left({t}\right)}\beta{\left({t}\right)}{B}{\left({t}\right)}={1},$

where $\kappa$ ,   $\theta$ , $\sigma$ , $\alpha$ and   $\beta$ are deterministic known functions. The system is subject to terminal conditions   ${A}{\left({T}\right)}={0}$ and ${B}{\left({T}\right)}={0}.$

In wasora, this problem can be solved by changing the terminal conditions to initial conditions with a change of variable $\tau={T}-{t}$ .

In that context, the Riccati system becomes:

$\frac{{{d}{A}{\left(\tau\right)}}}{{{\left.{d}{t}\right.}}}-\kappa{\left({T}-\tau\right)}\theta{\left({T}-\tau\right)}{B}{\left(\tau\right)}+\frac{1}{2}{\sigma}^{2}{\left({T}-\tau\right)}\alpha{\left({T}-\tau\right)}{B}{\left(\tau\right)}={0},$

$-\frac{{{d}{B}{\left(\tau\right)}}}{{{\left.{d}{t}\right.}}}-\kappa{\left({T}-\tau\right)}{B}{\left(\tau\right)}+\frac{1}{2}{\sigma}^{2}{\left({T}-\tau\right)}\beta{\left({T}-\tau\right)}{B}{\left(\tau\right)}={1},$

subject to initial conditions:

${A}{\left({0}\right)}={0}$ and   ${B}{\left({0}\right)}={0}.$

In reality, I have to solve this system several times, for different values of   ${T}$ . This means that instead of obtaining only one solution for   ${A}{\left(\tau\right)}$ and   ${B}{\left(\tau\right)}$ (which can be transformed to   ${A}{\left({t}\right)}$ and   ${B}{\left({t}\right)}$ ), I have to obtain several solutions   ${A}_{T}{\left(\tau\right)}$ and   ${B}_{T}{\left(\tau\right)}$ for each ${T}$ . Once I have them all, I will have   ${A}{\left(\tau,{T}\right)}$ and   ${B}{\left(\tau,{T}\right)}$ surfaces or, equivalently,   ${A}{\left({t},{T}\right)}$ and   ${B}{\left({t},{T}\right)}$ , which I will be able to interpolate.

Now that I have explained the problem, I would like to ask you how should I code this problem internally in wasora. Should I follow the steps that wasora performs when I use an input with these DAEs? Do you think there is a faster way to solve the system?

Of course, I am very familiar with wasora and once I follow the IDA steps it performs I will be able to code this properly, but I wanted to check with you for better ideas.

Thanks and have a nice weekend!

Rami

Jeremy Theler

unread,

Apr 24, 2019, 2:58:53 PM4/24/19

to was...@seamplex.com

So this is about speed, not about the possibility to solve the equations with wasora, right?

First, it seems that this equation can be solved with wasora using an input file (i.e. no external coding).

If you do that, can you add your input file to the examples directory?

Then, about speed. If you really need to solve this many times as fast as possible, do not guess but profile before optimizing.

Use valgrind and kcachegrind to profile. My guess (as I think it its yours) is that the bottle neck would be the evaluation of the residual.

If you confirm that, replace the general evaluation of the residual with a hard-coded version of the ODEs.

I think that is the fastsest thing you can do.

Second, instead of using IDA try KINSOL. The first one is for DAEs and the second one is for ODEs. Being more particular, maybe it would be faster too.

Third, I only tested wasora with SUNDIALS v2.x and now they are in v4.x. I think the API changed, so I would suggest you try to compile wasora against the newest version (and push the fixeds!)

--

jeremy

--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+un...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.
To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/35313D84-72DD-4E28-AD5B-EE70C9A0C1E5%40gmail.com.
For more options, visit https://groups.google.com/a/seamplex.com/d/optout.

Jeremy Theler

unread,

May 1, 2019, 5:45:42 AM5/1/19

to was...@seamplex.com

Ramiro, when you have time take a look at the new SUNDIALS 3.x

They allow for choosing the linear solver to approximate the nonlinear solution.

Maybe tweaking that kind of stuff will lead to faster solutions.

There are also SUNDIALS 4.x and even an alpha 5.x

When you do take a look, merge your code into wasora!

--

jeremy

Ramiro Vignolo

unread,

May 5, 2019, 4:52:37 PM5/5/19

to was...@seamplex.com

Sorry for the delay.

Regarding your first email:

1. I have already uploaded a Riccati ODEs example to wasora.

2. The evaluation of the residuals would be the bottleneck. I have hardcoded those residuals but they still depend on some expression that explicitly depend on time, so wasora_evaluate_expression will be called. However, I do not think this will be a problem.

3. Did you mean KINSOL or CVODES?

Regarding your second email:

1. Okay, I will make time to take a look at new SUNDIALS versions.

2. I have some comments with respect to the current implementation of transient loops when DAEs are present. Please, take a look at this issue.

Thanks,

Rami

On 1May, 2019, at 06:45, Jeremy Theler <jer...@seamplex.com> wrote:

Ramiro, when you have time take a look at the new SUNDIALS 3.x
They allow for choosing the linear solver to approximate the nonlinear solution.
Maybe tweaking that kind of stuff will lead to faster solutions.

There are also SUNDIALS 4.x and even an alpha 5.x

When you do take a look, merge your code into wasora!

--
jeremy

On Wed, 2019-04-24 at 15:58 -0300, Jeremy Theler wrote:

So this is about speed, not about the possibility to solve the equations with wasora, right?
First, it seems that this equation can be solved with wasora using an input file (i.e. no external coding).
If you do that, can you add your input file to the examples directory?

Then, about speed. If you really need to solve this many times as fast as possible, do not guess but profile before optimizing.
Use valgrind and kcachegrind to profile. My guess (as I think it its yours) is that the bottle neck would be the evaluation of the residual.
If you confirm that, replace the general evaluation of the residual with a hard-coded version of the ODEs.
I think that is the fastsest thing you can do.

Second, instead of using IDA try KINSOL. The first one is for DAEs and the second one is for ODEs. Being more particular, maybe it would be faster too.

Third, I only tested wasora with SUNDIALS v2.x and now they are in v4.x. I think the API changed, so I would suggest you try to compile wasora against the newest version (and push the fixeds!)

--
jeremy

On Fri, 2019-04-19 at 19:38 -0300, Ramiro Vignolo wrote:

Hi all! It's been a while...

As you may recall, sometimes I come around and ask for your help in some kind of, what I think, is a nice problem. The last time was about a puzzle (which I have solved XD).

Let's focus on today's problem. I have to solve the following system of Riccati ODEs:

<send-math-email_6162.png>

<send-math-email_6163.png>

where <send-math-email_6183.png>, <send-math-email_6128.png>, <send-math-email_6184.png>, <send-math-email_6185.png> and <send-math-email_6165.png> are deterministic known functions. The system is subject to terminal conditions <send-math-email_6160.png> and <send-math-email_6166.png>

In wasora, this problem can be solved by changing the terminal conditions to initial conditions with a change of variable <send-math-email_6167.png>.

In that context, the Riccati system becomes:

<send-math-email_6186.png>

<send-math-email_6169.png>

subject to initial conditions:

<send-math-email_6170.png> and <send-math-email_6171.png>

In reality, I have to solve this system several times, for different values of <send-math-email_6172.png>. This means that instead of obtaining only one solution for <send-math-email_6173.png> and <send-math-email_6174.png> (which can be transformed to <send-math-email_6175.png> and <send-math-email_6176.png>), I have to obtain several solutions <send-math-email_6177.png> and <send-math-email_6178.png> for each <send-math-email_6172.png>. Once I have them all, I will have <send-math-email_6179.png> and <send-math-email_6180.png> surfaces or, equivalently, <send-math-email_6181.png> and <send-math-email_6182.png>, which I will be able to interpolate.

Now that I have explained the problem, I would like to ask you how should I code this problem internally in wasora. Should I follow the steps that wasora performs when I use an input with these DAEs? Do you think there is a faster way to solve the system?

Of course, I am very familiar with wasora and once I follow the IDA steps it performs I will be able to code this properly, but I wanted to check with you for better ideas.

Thanks and have a nice weekend!

Rami

--
You received this message because you are subscribed to the Google Groups "wasora" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wasora+un...@seamplex.com.
To post to this group, send email to was...@seamplex.com.
Visit this group at https://groups.google.com/a/seamplex.com/group/wasora/.

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/9319ee5e9da12ab98f40e5949f34106817d9079c.camel%40seamplex.com.

Germán Theler

unread,

May 5, 2019, 6:13:22 PM5/5/19

to was...@seamplex.com

hardcode the function of t into the residual, avoid calling wasora_evaluate()

--
jeremy theler
www.seamplex.com

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/C55156F4-1B2C-44EF-B4F8-5CCF7D17C090%40gmail.com.

Ramiro Vignolo

unread,

May 5, 2019, 6:27:22 PM5/5/19

to was...@seamplex.com

I cannot do that because the residuals depend on expressions that must be provided by the user. Think of them as macroscopic cross sections in milonga. One of the coolest thing about wasora/milonga is that those cross sections can be defined as functions (algebraic, point wise, etc) or, more generally, algebraic expressions.

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/CAK0Liy%3DrpfTvF0Q-ctK%3DaxgcT7LrmJV08PRyysd1qkPPrT9f5A%40mail.gmail.com.

Germán Theler

unread,

May 5, 2019, 6:29:57 PM5/5/19

to was...@seamplex.com

then ask the user to define functions with a certain name, then acces the function_t internals from your residual

you can acces the pointwise data there

--
jeremy theler
www.seamplex.com

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/40ED81CF-BA3B-47D8-87C6-7AA0527F2478%40gmail.com.

Ramiro Vignolo

unread,

May 5, 2019, 7:00:33 PM5/5/19

to was...@seamplex.com

You want me to use wasora_evaluate_function() instead of wasora_evaluate_expression() in the computation of residuals, right?

The problem is that I do not want to use user defined functions only, but I want to be able to use user defined expressions. Expressions are more general and they can be made up using user defined functions and everything else.

The approach that you are pointing out is the one that we used for pcex back in the old days, right? However, this is not what I want, this is only one possible case. Let me show you the instruction I am coding:

SHORT_RATE NAME r1 TYPE affine {

KAPPA kappa_r1(t)

THETA theta_r1(t)

SIGMA sigma_r1(t)

ALPHA alpha_r1(t)

BETA beta_r1(t)

...

}

This instruction will solve several ODEs system for different terminal (or initial) conditions and built surfaces A(t,T) and B(t,T) that are needed for posterior instructions. As you see, kappa_r1(t) should be a function previously defined. For example:

FUNCTION kappa_r1(t) FILE_PATH kappa_r1.dat INTERPOLATION linear

However, in a more general case, I would like to be able to handle the following situation:

SHORT_RATE NAME r1 TYPE affine {

KAPPA (1 / (kappa_r1_b - kappa_r1_a)) * integral(kappa_r1(t), t, kappa_r1_a, kappa_r1_b)

THETA theta_r1(t)

SIGMA sigma_r1(t)

ALPHA alpha_r1(t)

BETA beta_r1(t)

...

}

At the moment, I am following your previous recommendation: "do not guess but profile before optimizing."

Thanks for your insight!

Rami

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/CAK0LiykBn4hfbZp9w%3DUMGF5H%3D4%2BTXsb8uPbpqbofYON_y7xgPQ%40mail.gmail.com.

Germán Theler

unread,

May 5, 2019, 7:26:01 PM5/5/19

to was...@seamplex.com

efficiency comes at the price of complexity

do this:

1. in a first input allow for general expressions and dump the data for the time range of interest into a file

2. in a second input read back the funci

tions as pointwise data

in the evaluation of the residuals keep a static index of the last interval for which the independent variable of the data contain the current time, and use it to interplate linearly. When time advances enough, add one to the index (in some rare cases time might go backward if the residual is big and dt needs to be reduced, take this into account)

this scheme is faster than having to bisect each time step to find the pair of datums to interpolate

--
jeremy theler
www.seamplex.com

To view this discussion on the web visit https://groups.google.com/a/seamplex.com/d/msgid/wasora/7441EE60-EA55-4C4A-AAFA-FB2D1297C139%40gmail.com.

jeremy theler

unread,

May 5, 2019, 7:30:19 PM5/5/19

to was...@seamplex.com

no need to use two runs and files you can do 1 during the static step and write the pointiwse data into internal arrays and then interpolate that cached information as in 2 in each step

--
jeremy theler
www.seamplex.com

Reply all

Reply to author

Forward