EME161 -- Combustion and the Environment

#### A sample program in C/C++ interpreter Ch

/********************************************************************* * 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(); }

#### A sample program in MATLAB

% 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

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