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,

Please find the code below:
#include
#include
#include
#include "intel_ode.h"
#include "mkl.h"
using namespace std;

double resistance=-40.8;
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
double mat[25];

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

mat[0]=-c;
mat[1]=c;
mat[2]=0.0;
mat[3]=0.0;
mat[4]=0.0;

mat[5]=d;
mat[6]=-2.172-d;
mat[7]=2.172;
mat[8]=0.0;
mat[9]=0.0;

mat[10]=0.0;
mat[11]=1.077;
mat[12]=-2*a-1.077;
mat[13]=a;
mat[14]=a;

mat[15]=0.0;
mat[16]=0.0;
mat[17]=b;
mat[18]=-b-f;
mat[19]=f;

mat[20]=0.0;
mat[21]=0.0;
mat[22]=g;
mat[23]=e...

from ia32 folder for my 32-bit laptop).

Regards,
-Sergey G

[ted_kord]

Hi Sergey

I'll add it to my toolbox. Thank you.

Ted

[Sergey_G_Intel]

Good luck, Ted!

-Sergey G