Specific Stiff ODE Problem

[ted_kord]

Hi all

I have a set of five coupled ODEs.
Each ode has one or more additional parameters.
Each parameter however, is a function of ,say, voltage.
So, for example:
ODE 1: RHS = f(t,y) = a * exp(......)
where the value of 'a' is dependent on a specific voltage, i.e, a(voltage).

I need to (must) solve the ODEs and obtain the current solution at each time step!

Question 1 : How can I supply at each time step a different voltage to calculate the parameter 'a'?
Question 2a: How can I obtain the solution at each time step ?
Question 2b: Perform some operation(s) with the current solution at that time step and then advance to the next time step?

I've tried several stiff ODE algorithms but no luck. I'd now like to use the Intel ODE library to solve this problem.

Any pointers/help in how to do this is appreciated.

Thanks

Ted Kord

[Sergey_G_Intel]

Hi Ted!

First of all, thank you for trying Intel ODE Solver Library! To answer your questions, let me consider general Fortran interface CALL dodesol(ipar,n,t,t_end,y,rhs,jacmat,h,hm,ep,tr,dpar,kd,ierr).

According to the Manual, you can set parameter ipar(3) to 1. This would slow down the solver a bit, but it will exit after each step. As a result, you will be able to do your operations and/or do whatever you like to the computed solution after each step. You can find current time in t and current solution in y. This could be summarized as follows:
ipar(3)=1
do while (tCALL dodesol(ipar,n,t,t_end,y,rhs,jacmat,h,hm,ep,tr,dpar,kd,ierr)
! do what you need to do with the computed solution y at time t
enddo
This should answer questions 2a and 2b.

I dont fully get your question number 1. Do you mean that you would like to have a way to pass additional parameter voltage to RHS routine? If so, then I would need more details about your problem to help you. I can think of stuff like global variables/common blocks or routine like get_current_voltage() to get the information inside RHS routine without passing it through fixed set of parameters.

Regards,
Sergey G

[ted_kord]

Hi Sergey,

Thanks for the help. Your explanation makes sense 'BUT':

I've set ipar[3] = 1.
I carried out the calcluations I required after the dodesol* call.

However, it seems to be taking an inordinately long time before completing a step (somewhat understably).
I actually have to stop the run because the time taken to compute just that step is totally unacceptable for the purposes of my simulation.

Am I doing something wrong?

Example code:

...
...
ipar[3] = 1

do {
dodesol(ipar, &n, &t, &t_end, y, rhs_v_d_p, jacmat_v_d_p, &h, &hm, &ep, &tr, dpar, kd, &ierr);

//Perform some calcultions
.......
.......

//Output results
cout ....................

} while (t < t_end);

[Sergey_G_Intel]

Hi Ted,

What I can see from your explanations is that your are using C code. C arrays have slightly different indexing than Fortran ones (the latter indexing should be shifted by -1 to get the former one). Therefore, you have to set ipar[2]=1 to exit after each step of computations. With ipar[3] set to 1, the dodesol routine would try to compute your Jacobi matrix using provided by you routine jacmat_v_d_p and compute the solution at time t_end. This may take a while if your problem is stiff.
In summary, try to set ipar[2] to 1 and if this does not help much, then I would need more details about your test case.

Regards,
Sergey G

[ted_kord]

You're absolutely spot on Sergey. I just remembered that now.

It's working fine now. Thx.


I'll keep you updated.

Ted

[Sergey_G_Intel]

You're always welcome, Ted!

-Sergey G

[ted_kord]

Hi Sergey

The program runs fine, performs the intermediate calculatations and outputs fine at the required steps. 'BUT' now I have a different question/problem.

The values output are extremely large which shouldn't be the case. However, I get somewhat similar values with one or two other stiff integrators. ~ 1e+305

Then after a certain time step, (less than 1% of the total time interval) , it spits out an 'inf' which cnonsequently renders other calculations as 'nan'. Then, the program just hangs.

Any thoughts?

Ted

[Sergey_G_Intel]

Hi Ted,

Numbers like 1e+305 usually appears as a result of usage of uninitialized variables that contain garbage or memory corruption, e.g., reading u[100] in array u of size 100. Can I get a test case (preferably, small) that demonstrates your problem?

Regards,
-Sergey G

[ted_kord]

Hi Sergey

Here's a test case. This is the set of ODEs I'm trying to solve:

//Constants
a = 36.367689741213006 * exp(0.055398582982867296 * resistance);
b = 0.28447798136999997 * exp(-0.02159796106832056 * resistance);
c = 176.39278458673601 * exp(0.055398582982867296 * resistance);
d = -4.6565444211489995 * exp(-0.036607982538678523 * resistance);
e = 1334.3934962347917 * exp(-0.030581515065542283 * resistance);
f = 5003.889182919519 * exp(0.018865463394165619 * resistance);
g = e*b/f;

//ODEs
dY1 = (d * Y2)-(c * Y1)
dY2 = -((d + 2.172) * Y2)+(c * Y1)+(1.077 * Y3)
dY3 = -((1.077 + a + a) * Y3) + (2.172 * Y2) + (b * Y4) + (g * Y5)
dY4 = -((b + f) * Y4) + (a * Y3) + (e * Y5)
dY5 = -((g + e) * Y5) + (a * Y3) + (f * Y4)

Ideally, 'start_time' = 0.0 and 'end_time' = 6100.0 and over. 'Time_step' = 0.02 or less.
Again, ideally, I'd like to obtain solutions at each time step.

N.B: Resistances have values of -80 to 40 in steps of 10. However, a test
resistance of -40 should be sufficient.

Thanks,

Ted

[Sergey_G_Intel]

Hi Ted,

Do you have initial condition Y(t=0)?

Regards,
-Sergey G

[ted_kord]

Y1 = Y2 = Y3 = Y4 = Y5 = 0.2;

Although, this is not set in stone.

Arbitrary values less than 1 are okay.

Ted

[Sergey_G_Intel]

Hi Ted,

Ive created a small code that computes the solution to your ODE system:

#include
#include
#include
#include "intel_ode.h"
using namespace std;

double resistance=-40.0;
double a = 36.367689741213006 * exp(0.055398582982867296 * resistance);
double b = 0.28447798136999997 * exp(-0.02159796106832056 * resistance);
double c = 176.39278458673601 * exp(0.055398582982867296 * resistance);
double d = -4.6565444211489995 * exp(-0.036607982538678523 * resistance);
double e = 1334.3934962347917 * exp(-0.030581515065542283 * resistance);
double f = 5003.889182919519 * exp(0.018865463394165619 * resistance);
double g = e*b/f;

void rhs (int *n, double *t, double *y, double *rhs) {
rhs[0] = (d * y[1])-(c * y[0]);
rhs[1] = -((d + 2.172) * y[1])+(c * y[0])+(1.077 * y[2]);
rhs[2] = -((1.077 + a + a) * y[2]) + (2.172 * y[1]) + (b * y[3]) + (g * y[4]);
rhs[3] = -((b + f) * y[3]) + (a * y[2]) + (e * y[4]);
rhs[4] = -((g + e) * y[4]) + (a * y[2]) + (f * y[3]);
}

int main(){

int ipar[128];
int n=5, ierr;
int kd[5];
double t=0.0, t_end=6100.0, ti=0.02;
double y[5], h=1.0e-7, hm=1.0e-14, ep=1.0e-3, tr=1.0e-3, *dpar;
void *jacmat=NULL;//dummy at the moment

cout<<"Starting computations...n";
for (int i=0;i<128;i++) ipar=0;

dpar=(double*)malloc((7+2*n)*n*sizeof(double));
//ipar[2]=1; //Commented out as solver exits for t=0.02*k anyway

for (int i=0;i<5;i++) y=0.2;

int steps=(int)(t_end/ti);
cout<<"Performing "<<<" steps with dodesoln";
for (int k=0;kt=k*ti;
t_end=(k+1)*ti;
do {
dodesol(ipar,&n,&t,&t_end,y,rhs,jacmat,&h,&hm,&ep,&tr,dpar,kd,&ierr);
if (ierr) {
cout<<"dodesol exited with error "<<<"n";
}
} while (tif (!((k+1)%5000)) {
cout<<"At time t="<<<" the following solution obtained:n";
cout<<"y[0]="< cout<<"y[1]="< cout<<"y[2]="< cout<<"y[3]="< cout<<"y[4]="<}
}

system("PAUSE");
return 0;
}

It computes the solution at time tk=0.02*k, k=1,,305000. It prints the solution at time tj=100.0*j, j=1,..,61. In a couple of seconds on my laptop, Im getting the solution at time t=6100.0.
At time t=6100 the following solution obtained:
y[0]=-0.052425
y[1]=0.050077
y[2]=0.100991
y[3]=0.593447
y[4]=0.30791

As you can see, there are no NaNs or numbers like 1.0e+305. Could you please verify that the program computes the solution you need? Thank you!

Regards,
-Sergey G

[ted_kord]

Hi Sergey

Thanks for this. Works fine at resistance = -40.0. However, when the resistance is set to -80, it hangs. Same thing happens at other resistance values, i.e. -70, etc.

I'm tearing my hair out:)

Thanks

Ted

[ted_kord]

Hi Sergey

Update.
It seems like any resistance vlaues more positive than -40, i.e -39, -38 , .... work fine. Resistance more negative than -40, i.e -41, -42, ....-50,....... don't.

Thanks,

Ted

[Sergey_G_Intel]

Hi Ted,

It looks like your system has a solution growing up to infinity for resistances <=-41. You mentioned once that your real system has variable parameters. May be, this is the only way to keep solution from growing too far, e.g., having resistance <-40 for a short period of time. Besides, I wonder why the equation for Y3 has the form:
dY3 = -((1.077 + a + a) * Y3) + (2.172 * Y2) + (b * Y4) + (g * Y5)
with (..+a+a). Anything special about it?

Anyway, I'll try to lookif anythingcan be done here. You can verify once again that the system is correct and that it can stay with resistances below -40 for a long while.

Regards,
-Sergey G

[ted_kord]

Hi Sergey

Thanks. If anything can be done for the case, resistance < -40.8, I shall be eternally in your debt. I've been working at this for two straight weeks and it finally seems like there's light at the end of the tunnel.

No, there's nothing special about the equation for Y3. The (..a+a) arose during factorisation and I just decided to leave it like that (for no special reason). (2*a) is much more concise anyway.

Meanwhile, I'll keep on running tests for resistance >= -40

Thanks

Ted

[Sergey_G_Intel]

Hi Ted,

It looks to me that your system is getting a positive eigenvalue (or an eigenvalue with positive real part, to be precise) around resistance -40.9. As the system is relatively simple, I'll try to check this in the near future. If this is the case, then there is no way for you to get bounded solution, sorry...

Regards,
-Sergey G

[Sergey_G_Intel · Accepted solution]

Hi Ted,

Ive estimated eigenvalues of your problem using Intel MKL LAPACK. As a result, for resistance=-40.8 Ive got that the eigenvalues of your system are:

0: -9.32634+0i
1: -0.0462603+6.01633i
2: -0.0462603+-6.01633i
3: -1.06398e-013+0i
4: -6965.59+0i

As you can see, the real parts of all eigenvalues are negative. This means that the system has bounded solution. When I set resistance to -40.9, Ive got the following set of eigenvalues:

0: -9.2891+0i
1: 0.0431913+5.99626i
2: 0.0431913+-5.99626i
3: 5.6765e-013+0i
4: -6975.47+0i

As you can see, there are 2 eigenvalues with positive real parts. Asymptotically, the solution of the system would behave as exp(0.043*t), which is unbounded function when t grows. This is what you see when trying to solve the system with resistance <-40.8 with our ODE solver.

It looks that I cannot help you more with your ODE system as it is just unstable for certainvalues of parameter "resistance". Sorry, if this does not help you much

Regards,
-Sergey G

[ted_kord]

Hi Sergey

That's very useful information. Thanks anyway for all the help. It's been very enlightening.

Keep up the good work.

Best regards

Ted

[Sergey_G_Intel]

Your are always welcome, Ted!

Good luck!

-Sergey G

[ted_kord]

Hi Sergey

Thanks again for all your help earlier.

Could you please provide a copy of the code you used to estimate the eigenvalues of my problem using Intel MKL LAPACK. I anticipate needing to do this quite a few times in the future.

Thanks

Ted