EME161 -- Combustion and the Environment
A sample program in C/C++ interpreter Ch
<![CDATA[
/*********************************************************************
* File: command.ch *
* This C/C++ program solves an ODE with multiple variables and plot *
* the results. To run this program, use computers in CAE lab or *
* download * C/C++ interpreter Ch Professional Edition from *
* http://www.softintegration.com/download *
* ******************************************************************/
#include <numeric.h>
#include <chplot.h>
#define NVAR 2 // number of variables
#define POINTS 256 // number of points
// Global variables passed to the ode function
double p =1, r = 1, eps = 0.0001;
int heaviside(double x) {
if(x < 0)
return 0;
else
return 1;
}
void func(double t, double y[], double dydt[], void *param) {
dydt[0] = -p*heaviside(y[0]-eps)*y[1]/(pow(y[0], 0.3)); // ode 1
dydt[1]= -r*y[0]*y[1]; // ode 2
}
int main() {
// Initial condition vector
double y0[NVAR] = {1, 1};
// Total integration time interval
double t0=0, tf=2;
// variables
double t[POINTS], y[NVAR][POINTS];
// plotting class
class CPlot plot;
// use Runge-Kutta method to solve the ODE
oderk(func, t0, tf, y0, NULL, t, y);
// plotting
plot.data2D(t, y);
plot.title("The solution for the ODE");
plot.label(PLOT_AXIS_X, "time (seconds)");
plot.label(PLOT_AXIS_Y, "y1 and y2");
plot.legend("y1", 0);
plot.legend("y1", 1);
plot.plotting();
}
]]>
Output
A sample program in MATLAB
<![CDATA[
% Place the following text in the Matlab command window
y0 = [1 1];
% Initial condition vector
tspan = [0 2];
% Total integration time interval
global P R EPS
% Variables to pass to the ode function
P = 1; R = 1;
% Set parameter values
EPS = 0.0001;
% Set parameter values
[t,y] = ode15s(@odetest, tspan, y0);
% Find the solution
plot(t,y(:,1),'-',t,y(:,2),'-')
% Plot the solution
legend('y1','y2');
% Plot legend
xlabel('time');
% Plot x-axis label
ylabel('y1, y2');
% Plot y-axis label
]]>
A sample program in MATLAB
<![CDATA[
function dydt = odetest(t,y)
% Defines the ode equation set function
global P R EPS
% Variables passed to the ode equation set
dydt = zeros(2,1);
% Defines the ode vector
dydt(1) = -P*heaviside(y(1) - EPS)*y(2)/(y(1)^0.3); % ode 1
dydt(2) = -R*y(1)*y(2); % ode 2
]]>