Results 1 to 3 of 3

Thread: Differential Equations Simulation

  1. #1
    Join Date
    Jul 2012
    Posts
    1
    Qt products
    Qt4

    Default Differential Equations Simulation

    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:

    Qt Code:
    1. #include <iostream>
    2. #include "nr3.h"
    3. #include "odeint.h"
    4. #include "ludcmp.h"
    5. #include "stepper.h"
    6. #include "stepperross.h"
    7. #include <fstream>
    8. #include <cmath>
    9. using namespace std;
    10.  
    11.  
    12. double a = 0.3646; // parameters for equation (1) in paper
    13. double b = 8.6971;
    14. double k1 = 2.0732;
    15. double c = 0.4467;
    16.  
    17. double alpha = 10.0; // parameters for equation (2) in paper
    18. double beta = 10.0;
    19. double k2 = 0.2041;
    20. double S = 0.1;
    21.  
    22. double gamma = 0.77; // "rise factor" (nM)
    23. double eta = 0.91; // "decay rate" (/s)
    24. double t_ox = 60.0; // "time constant for orexin" (s)
    25. double t_ser = 10.0; // "time constant for serotonin/5_HT" (s)
    26.  
    27. double phi = 33.57; //"Release per stimulus rate" (nM)
    28. double Km = 170.0; // (nM) Km Vmax Michaelis-Menton parameters
    29. double Vmax = 1300; // (nM/s)
    30.  
    31. //NB Input/output (3) and (4) were substiuted directly subsitituted into the equations (1) and (2) in the paper
    32. //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))) ;
    33. // = A1*y[0] + A2 + A3*1/(1+ A4*pow(y[1],A5))
    34.  
    35. double A1 = (-1/t_ox);
    36. double A2 = (a/t_ox);
    37. double A3 = (b/t_ox);
    38. double A4 = exp(k1/c);
    39. double A5 = (-1/c*log(10.0));
    40.  
    41.  
    42. //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)) ;
    43. // = B1*y[2] + B2 + B3* 1/(1+B4*pow(y[3],B5))
    44.  
    45. double B1 = (-1/t_ser);
    46. double B2 = (alpha/t_ser);
    47. double B3 = (beta/t_ser);
    48. double B4 = exp(k2/S);
    49. double B5 = (-1/S*log(10.0));
    50.  
    51. //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);
    52. // = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2);
    53.  
    54. double C1 = (b*A4/c*log(10.0)*t_ox);
    55. double C2 = A5-1;
    56. //double C3 = 1 + A4*pow(y[1],A5);
    57.  
    58.  
    59. //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)
    60. // = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2);
    61.  
    62. double D1 = (beta*B4/S*log(10.0)*t_ser);
    63. double D2 = B5-1;
    64. //double D3 = 1 + B4*pow(y[3],B5);
    65.  
    66. //==========================================================================================================================
    67.  
    68. struct rhs
    69. {
    70.  
    71. void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
    72. {
    73. dydx[0] = A1*y[0] + A2 + A3*1/(1 + A4*pow(y[1],A5)) ;
    74. dydx[1] = -eta*y[1] + gamma*y[2];
    75. dydx[2] = B1*y[2] + B2 + B3*1/(1 + B4*pow(y[3],B5)) ;
    76. dydx[3] = phi*y[0] - Vmax*y[3]/(Km + y[3]);
    77. }
    78.  
    79. void jacobian(const Doub x, VecDoub_I &y, VecDoub_O &dfdx,
    80. MatDoub_O &dfdy)
    81. {
    82. Int n=y.size();
    83. for (Int i=0;i<n;i++)
    84.  
    85. dfdx[i]=0.0;
    86. dfdy[0][0] = A1;
    87. dfdy[0][1] = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2));
    88. dfdy[0][2] = 0.0;
    89. dfdy[0][3] = 0.0;
    90. dfdy[1][0] = 0.0;
    91. dfdy[1][1] = -eta;
    92. dfdy[1][2] = gamma;
    93. dfdy[1][3] = 0.0;
    94. dfdy[2][0] = 0.0;
    95. dfdy[2][1] = 0.0;
    96. dfdy[2][2] = B1;
    97. dfdy[2][3] = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2));
    98. dfdy[3][0] = phi;
    99. dfdy[3][1] = 0.0;
    100. dfdy[3][2] = 0.0;
    101. dfdy[3][3] = -Vmax*Km/pow(Km+y[3],2);
    102.  
    103. }
    104. };
    105.  
    106.  
    107.  
    108. int main()
    109. {
    110. ofstream outFile("data.txt"); // Opens and instantiate a data output file
    111.  
    112. const Int nvar = 4; //Number of equations to compute
    113. 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)
    114. VecDoub ystart(nvar);
    115.  
    116. //Initial conditions
    117. ystart[0]=0.5; // f_DRN(t=0) (Hz)
    118. ystart[1]=2.8; // [Ox](t=0) (nM)
    119. ystart[2]=5.0; // f_LHA(t=0) (Hz)
    120. ystart[3]=1.6; // [Ser](t=0) (nM)
    121.  
    122. Output out(500); //the value inside out is the user defined number of points required!
    123. rhs d;
    124.  
    125. Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
    126. ode.integrate();
    127.  
    128. for (Int i=0;i<out.count;i++)
    129. {
    130. //output to console window
    131. cout << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
    132. out.ysave[2][i] << " " << out.ysave[3][i] << endl;
    133.  
    134. //write data to data.txt
    135. outFile << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
    136. out.ysave[2][i] << " " << out.ysave[3][i] << endl;
    137.  
    138. }
    139.  
    140. cout << endl << out.count << endl; //displays the number of points output by the program
    141.  
    142. outFile.close(); // Closes the data output file!
    143.  
    144. system("Pause");
    145.  
    146. }
    To copy to clipboard, switch view to plain text mode 

  2. #2
    Join Date
    Sep 2011
    Location
    Manchester
    Posts
    538
    Thanks
    3
    Thanked 106 Times in 103 Posts
    Qt products
    Qt4 Qt/Embedded
    Platforms
    MacOS X Unix/X11 Windows

    Default Re: Differential Equations Simulation

    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

  3. #3
    Join Date
    Feb 2006
    Location
    Munich, Germany
    Posts
    3,278
    Thanked 872 Times in 820 Posts
    Qt products
    Qt3 Qt4 Qt/Embedded
    Platforms
    MacOS X Unix/X11 Windows

    Default Re: Differential Equations Simulation

    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

Similar Threads

  1. Best way to design a simulation class?
    By Yes in forum Qt Programming
    Replies: 10
    Last Post: 28th November 2011, 00:19
  2. Replies: 2
    Last Post: 8th February 2011, 14:15
  3. alt enter simulation in qstring
    By Dilshad in forum Newbie
    Replies: 6
    Last Post: 29th December 2010, 05:59
  4. physics simulation
    By scrasun in forum General Programming
    Replies: 6
    Last Post: 11th June 2009, 02:19
  5. How to solve simultanious equations...??
    By joseph in forum General Programming
    Replies: 2
    Last Post: 28th May 2008, 10:53

Bookmarks

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •  
Digia, Qt and their respective logos are trademarks of Digia Plc in Finland and/or other countries worldwide.