Differential Equations Simulation

[GGH84]

Hi, I'm working on a GUI for a system of differential equations. The idea is to create a GUI similar to the "Oscilloscope" in the Qwt example folder that will allow the user to tweak various parameters of the system and view the realtime affects of the simulation. I've already created a basic console and file output implementation using C++ NUMERICAL RECIPES routines, but I'm not sure how to go about implementing a GUI representation of the basic console application. Any advise you guys could give me to speed up my progress on this project would be greatly appreciated.

Here's the basic code made for the Numerical recipes routine:

#include <iostream>
#include "nr3.h"
#include "odeint.h"
#include "ludcmp.h"
#include "stepper.h"
#include "stepperross.h"
#include <fstream>
#include <cmath>
using namespace std;
 
 
double a = 0.3646;     //  parameters for equation (1) in paper
double b = 8.6971;
double k1 = 2.0732;
double c = 0.4467;
 
double alpha = 10.0;   //  parameters for equation (2) in paper
double beta = 10.0;
double k2 = 0.2041;
double S = 0.1;
 
double gamma = 0.77;   // "rise factor" (nM)
double eta = 0.91;	   // "decay rate" (/s) 
double t_ox = 60.0;	   // "time constant for orexin" (s)
double t_ser = 10.0;   // "time constant for serotonin/5_HT"  (s) 
 
double phi = 33.57;      //"Release per stimulus rate" (nM)
double Km = 170.0;		// (nM)   Km Vmax Michaelis-Menton parameters
double Vmax = 1300;     // (nM/s)
 
//NB  Input/output  (3) and (4) were substiuted directly subsitituted into the equations (1) and (2) in the paper
//dydx[0] = (-1/t_ox)*y[0] + (a/t_ox) + (b/t_ox)*(1/(1 + exp(k1/c)*pow(y[1],-1/c*log(10.0))) ; 
//        = A1*y[0] + A2 + A3*1/(1+ A4*pow(y[1],A5))
 
double A1 = (-1/t_ox);
double A2 = (a/t_ox);
double A3 = (b/t_ox);
double A4 = exp(k1/c);
double A5 = (-1/c*log(10.0));
 
 
//dydx[2] = (-1/t_ser)*y[2] + (alpha/t_ser) + (beta/t_ser)* 1/(1 + exp(k2/S)*pow(y[3],-1/S*log(10.0)) ;
//		  = B1*y[2] + B2 + B3* 1/(1+B4*pow(y[3],B5))	
 
double B1 = (-1/t_ser);
double B2 = (alpha/t_ser);
double B3 = (beta/t_ser);
double B4 = exp(k2/S);
double B5 = (-1/S*log(10.0));
 
//dfdy[0][1] = (b*exp(k1/c)/c*log(10.0)*t_ox)*pow(y[1],-1/c*log(10.0) - 1.0)/pow(1 + exp(k1/c)*pow(y[1], 1/c*log(10.0)),2);
//      	 = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2);
 
double C1 = (b*A4/c*log(10.0)*t_ox);
double C2 = A5-1;
//double C3 = 1 + A4*pow(y[1],A5);
 
 
//dfdy[2][3] = (beta*exp(k2/S)/S*log(10.0)*t_ser)*pow(y[3],-1/Slog(10.0) - 1)/pow(1+exp(k2/S)*pow(y[3],S*log(10.0)),2)
//           =  D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2);
 
double D1 = (beta*B4/S*log(10.0)*t_ser);
double D2 = B5-1;
//double D3 = 1 + B4*pow(y[3],B5);
 
//==========================================================================================================================
 
struct rhs 
{
 
	void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
	{
		dydx[0] = A1*y[0] + A2 + A3*1/(1 + A4*pow(y[1],A5)) ;
		dydx[1] = -eta*y[1] + gamma*y[2];
		dydx[2] = B1*y[2] + B2 + B3*1/(1 + B4*pow(y[3],B5)) ;
		dydx[3] = phi*y[0] - Vmax*y[3]/(Km + y[3]);
	}
 
void jacobian(const Doub x, VecDoub_I &y, VecDoub_O &dfdx,
MatDoub_O &dfdy)
{
	Int n=y.size();
	for (Int i=0;i<n;i++)
 
	dfdx[i]=0.0;
	dfdy[0][0] = A1;                       
	dfdy[0][1] = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2));
	dfdy[0][2] = 0.0;
	dfdy[0][3] = 0.0;
	dfdy[1][0] = 0.0;
	dfdy[1][1] = -eta;
	dfdy[1][2] = gamma;
	dfdy[1][3] = 0.0;
	dfdy[2][0] = 0.0;
	dfdy[2][1] = 0.0;
	dfdy[2][2] = B1;
	dfdy[2][3] = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2));
	dfdy[3][0] = phi;
	dfdy[3][1] = 0.0;
	dfdy[3][2] = 0.0;
	dfdy[3][3] = -Vmax*Km/pow(Km+y[3],2);
 
}
};
 
 
 
int main()
{
	ofstream outFile("data.txt");  // Opens and instantiate a data output file
 
	const Int nvar = 4;  //Number of equations to compute
	const Doub atol=1.0e-5,rtol=atol, h1=2.9e-4, hmin=0.0, x1=0.0, x2=250;  //temporal integration limits (x1->t1, x2->t2)
	VecDoub ystart(nvar);
 
	//Initial conditions
	ystart[0]=0.5;  // f_DRN(t=0) (Hz)
	ystart[1]=2.8;  // [Ox](t=0) (nM)
	ystart[2]=5.0;  // f_LHA(t=0) (Hz)
	ystart[3]=1.6;  // [Ser](t=0) (nM)
 
	Output out(500);  //the value inside out is the user defined number of points required!
	rhs d;
 
	Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
	ode.integrate();
 
for (Int i=0;i<out.count;i++)
{
	//output to console window
	cout << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " << 
	out.ysave[2][i] << " " << out.ysave[3][i] << endl;
 
	//write data to data.txt
    outFile << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " << 
	out.ysave[2][i] << " " << out.ysave[3][i] << endl;
 
}
 
cout << endl << out.count << endl;  //displays the number of points output by the program
 
	outFile.close();  //  Closes the data output file!
 
	system("Pause");
 
}

#include <iostream>
#include "nr3.h"
#include "odeint.h"
#include "ludcmp.h"
#include "stepper.h"
#include "stepperross.h"
#include <fstream>
#include <cmath>
using namespace std;


double a = 0.3646;     //  parameters for equation (1) in paper
double b = 8.6971;
double k1 = 2.0732;
double c = 0.4467;

double alpha = 10.0;   //  parameters for equation (2) in paper
double beta = 10.0;
double k2 = 0.2041;
double S = 0.1;

double gamma = 0.77;   // "rise factor" (nM)
double eta = 0.91;	   // "decay rate" (/s) 
double t_ox = 60.0;	   // "time constant for orexin" (s)
double t_ser = 10.0;   // "time constant for serotonin/5_HT"  (s) 

double phi = 33.57;      //"Release per stimulus rate" (nM)
double Km = 170.0;		// (nM)   Km Vmax Michaelis-Menton parameters
double Vmax = 1300;     // (nM/s)

//NB  Input/output  (3) and (4) were substiuted directly subsitituted into the equations (1) and (2) in the paper
//dydx[0] = (-1/t_ox)*y[0] + (a/t_ox) + (b/t_ox)*(1/(1 + exp(k1/c)*pow(y[1],-1/c*log(10.0))) ; 
//        = A1*y[0] + A2 + A3*1/(1+ A4*pow(y[1],A5))

double A1 = (-1/t_ox);
double A2 = (a/t_ox);
double A3 = (b/t_ox);
double A4 = exp(k1/c);
double A5 = (-1/c*log(10.0));


//dydx[2] = (-1/t_ser)*y[2] + (alpha/t_ser) + (beta/t_ser)* 1/(1 + exp(k2/S)*pow(y[3],-1/S*log(10.0)) ;
//		  = B1*y[2] + B2 + B3* 1/(1+B4*pow(y[3],B5))	

double B1 = (-1/t_ser);
double B2 = (alpha/t_ser);
double B3 = (beta/t_ser);
double B4 = exp(k2/S);
double B5 = (-1/S*log(10.0));

//dfdy[0][1] = (b*exp(k1/c)/c*log(10.0)*t_ox)*pow(y[1],-1/c*log(10.0) - 1.0)/pow(1 + exp(k1/c)*pow(y[1], 1/c*log(10.0)),2);
//      	 = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2);

double C1 = (b*A4/c*log(10.0)*t_ox);
double C2 = A5-1;
//double C3 = 1 + A4*pow(y[1],A5);


//dfdy[2][3] = (beta*exp(k2/S)/S*log(10.0)*t_ser)*pow(y[3],-1/Slog(10.0) - 1)/pow(1+exp(k2/S)*pow(y[3],S*log(10.0)),2)
//           =  D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2);

double D1 = (beta*B4/S*log(10.0)*t_ser);
double D2 = B5-1;
//double D3 = 1 + B4*pow(y[3],B5);

//==========================================================================================================================

struct rhs 
{

	void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
	{
		dydx[0] = A1*y[0] + A2 + A3*1/(1 + A4*pow(y[1],A5)) ;
		dydx[1] = -eta*y[1] + gamma*y[2];
		dydx[2] = B1*y[2] + B2 + B3*1/(1 + B4*pow(y[3],B5)) ;
		dydx[3] = phi*y[0] - Vmax*y[3]/(Km + y[3]);
	}

void jacobian(const Doub x, VecDoub_I &y, VecDoub_O &dfdx,
MatDoub_O &dfdy)
{
	Int n=y.size();
	for (Int i=0;i<n;i++)
		
	dfdx[i]=0.0;
	dfdy[0][0] = A1;                       
	dfdy[0][1] = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2));
	dfdy[0][2] = 0.0;
	dfdy[0][3] = 0.0;
	dfdy[1][0] = 0.0;
	dfdy[1][1] = -eta;
	dfdy[1][2] = gamma;
	dfdy[1][3] = 0.0;
	dfdy[2][0] = 0.0;
	dfdy[2][1] = 0.0;
	dfdy[2][2] = B1;
	dfdy[2][3] = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2));
	dfdy[3][0] = phi;
	dfdy[3][1] = 0.0;
	dfdy[3][2] = 0.0;
	dfdy[3][3] = -Vmax*Km/pow(Km+y[3],2);

}
};



int main()
{
	ofstream outFile("data.txt");  // Opens and instantiate a data output file

	const Int nvar = 4;  //Number of equations to compute
	const Doub atol=1.0e-5,rtol=atol, h1=2.9e-4, hmin=0.0, x1=0.0, x2=250;  //temporal integration limits (x1->t1, x2->t2)
	VecDoub ystart(nvar);
	
	//Initial conditions
	ystart[0]=0.5;  // f_DRN(t=0) (Hz)
	ystart[1]=2.8;  // [Ox](t=0) (nM)
	ystart[2]=5.0;  // f_LHA(t=0) (Hz)
	ystart[3]=1.6;  // [Ser](t=0) (nM)

	Output out(500);  //the value inside out is the user defined number of points required!
	rhs d;
	
	Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
	ode.integrate();

for (Int i=0;i<out.count;i++)
{
	//output to console window
	cout << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " << 
	out.ysave[2][i] << " " << out.ysave[3][i] << endl;

	//write data to data.txt
    outFile << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " << 
	out.ysave[2][i] << " " << out.ysave[3][i] << endl;

}

cout << endl << out.count << endl;  //displays the number of points output by the program

	outFile.close();  //  Closes the data output file!
	
	system("Pause");
	
}

To copy to clipboard, switch view to plain text mode 

[Spitfire]

What you want is fairly simple. It requires some nitty-gritty work but it's not a rocket science.

Build the UI you want, with qwt plot in the center and with bunch of spinboxes (or other controls) allowing to modify your parameters.
Connect to each control's 'change()' signal (there's specific change() signal emited by each control) to some 'recalculate' method which will re-do the math and replot the results at the end.
You're done

[Uwe]

Maybe it's worth to have a look at the FunctionData class of the sinusplot example ( sinusplot.cpp ). It shows a way how you can add the math only without having to build array of samples from it.

Uwe