- y''' + 5y'' + 8y' + 4y = 1

Use a timestep of dt = 200us and run the simulation for t = 10 seconds. The initial conditions are:

- y(0) = 5
- y'(0) = 1
- y''(0) = 0

This problem is fun because you can solve for a closed form solution and then compare that hand-solved answer to the simulated version. The simulated solution can be run using Matlab, C, or a combination. The combination solution is neat because it uses Matlab to create the data and plot the solution, but uses C-code (compiled in Matlab as a mex-function) as the super-fast solution engine. Mex-Functions are a bit tricky to learn but can often lead to valuable speed-ups in simulations.

The solution to the third order differential equation can be solved by hand. You should get the following function:

The plot for this function is shown below along with the iterated solution that I generated using C++:

**Simulated Solutions**

Here I present three numerical (i.e. simulated) solutions to the third order differential equation. All three yield the same simulated plot (see blue curve above).

**Matlab-Only Solution**

clear;clf;dt=200e-6;tMax=10;nSamples=tMax/dt;y=zeros(nSamples,1);a=1;b=0;y(1) =5;for i=2:nSamplesdy=a;da=b;db= -8*a-5*b-4*y(i-1) +1;y(i) =y(i-1) +dy*dt;a=a+da*dt;b=b+db*dt;endt= (1:length(y))*dt;plot(t,y);

Using the tic/toc commands in Matlab, I determined that the Matlab-only solution took an average of 12ms, including the time to create the plot.

**C++ Only Solution**

The C++ Only solution is fast, but we need a way to pass the data back to Matlab in order to plot it. In this solution, C++ writes the data to a binary file. Then a separate Matlab script reads the data from the file and plots it.

#include <iostream>

#include <cmath>

#include <fstream>

using namespace std;

int main()

{

double dt = 200e-6;

double tmax = 10;

double pi = 4 * atan(1);

int nSamples = floor(tmax/dt);

double da,db,dy;

double a=1;

double b=0;

double y[nSamples];

int i;

ofstream out("data2.bin",ios::out|ios::binary);

y[0] = 5;

for(i=1;i<nSamples;i++){

dy = a;

da = b;

db = -8*a - 5*b - 4*y[i-1] + 1;

y[i] = y[i-1] + dy*dt;

a = a + da*dt;

b = b + db*dt;

}

out.write((char *)&dt , sizeof(double));

out.write((char *)y ,nSamples*sizeof(double));

out.close();

return 0;

}

This is the Matlab plot code:

clear;clf

fid=fopen('data2.bin','rb');

dt=fread(fid,1,'double');y=fread(fid,'double');

fclose(fid);

t= (0:length(y)-1)*dt;

plot(t,y);xlabel('time (s)');title('Solution to 3rd Order Diff-Eq');

**C++ / Matlab / Mex Solution**

#include <mex.h>

#include <string.h>

#include <math.h>

// This is the subroutine that actually performs the simulation

void runSim(double **py, double **pt, double dt, double tMax, int *nSamples){

double pi = 4 * atan(1);

double da, db, dy;

double a=1;

double b=0;

int i;

double tTemp = dt;

double *y, *t;

*nSamples = floor(tMax/dt);

*py = new double[*nSamples];

y = *py;

*pt = new double[*nSamples];

t = *pt;

y[0] = 5;

t[0] = dt;

for(i=1;i<*nSamples;i++){

dy = a;

da = b;

db = -8*a - 5*b - 4* y[i-1] + 1;

y[i] = y[i-1] + dy*dt;

a = a + da*dt;

b = b + db*dt;

t[i] = t[i-1] + dt;

}

}

// ****************************************************

// ******************** START HERE ********************

// ****************************************************

// Mex routines must always start with "mexFunction"

// Here, the input data is imported from Matlab. Then the actual function

// executed (in this case "runSim"), and finally the output data is

// exported back to Matlab

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

// Step 1: Import input variables from Matlab

double dt = *(double *)mxGetData(prhs[0]);

double tMax = *(double *)mxGetData(prhs[1]);

// Step 2: Declare variables

int nSamples;

double *y,*t;

// Step 3: Run the simulation

runSim(&y, &t, dt, tMax, &nSamples);

// Step 4: Export the results back to Matlab

double *output;

if (nlhs>=1){

plhs[0] = mxCreateDoubleMatrix(nSamples, 1, mxREAL);

output = mxGetPr(plhs[0]);

memcpy(output, y, nSamples*sizeof(double));

}

if (nlhs>=2){

plhs[1] = mxCreateDoubleMatrix(nSamples, 1, mxREAL);

output = mxGetPr(plhs[1]);

memcpy(output, t, nSamples*sizeof(double));

}

// Step 5: Housekeeping

delete [] y;

delete [] t;

}

Once the C/C++ code has been written and compiled, the function is called from Matlab in the same way that any other function would be. In this case, the C file was called "diff_eq3.cpp", so therefore the compiled executable is called using "diff_eq3". This simulation took an average of 2ms, including plotting time. That is 6x faster than the Matlab-only solution!

clear; clf;

dt = 200e-6;

tmax = 10;

[y,t] = diff_eq3(dt,tmax);

plot(t,y);

Hi,

ReplyDeleteThank you so much for your great description of the three different kind of implementations side by side. I'm trying to create a mex-file and I started to learn from your code but when I was trying to compile it, I got some errors. It would be awesome if you can give some hints.

I'm using Matlab 7.12.0 (R2011), win 7 and Compiler: Microsoft Visual C++ 2010.

And here is what I got:

>> mex diff_eq3.c

diff_eq3.c

diff_eq3.c(14) : warning C4244: '=' : conversion from 'double' to 'int', possible loss of data

diff_eq3.c(15) : error C2065: 'new' : undeclared identifier

diff_eq3.c(15) : warning C4047: '=' : 'double *' differs in levels of indirection from 'int'

diff_eq3.c(15) : error C2143: syntax error : missing ';' before 'type'

diff_eq3.c(18) : error C2065: 'new' : undeclared identifier

diff_eq3.c(18) : warning C4047: '=' : 'double *' differs in levels of indirection from 'int'

diff_eq3.c(18) : error C2143: syntax error : missing ';' before 'type'

diff_eq3.c(54) : error C2143: syntax error : missing ';' before 'type'

diff_eq3.c(57) : error C2065: 'output' : undeclared identifier

diff_eq3.c(57) : warning C4047: '=' : 'int' differs in levels of indirection from 'double *'

diff_eq3.c(58) : error C2065: 'output' : undeclared identifier

diff_eq3.c(58) : warning C4022: 'memcpy' : pointer mismatch for actual parameter 1

diff_eq3.c(62) : error C2065: 'output' : undeclared identifier

diff_eq3.c(62) : warning C4047: '=' : 'int' differs in levels of indirection from 'double *'

diff_eq3.c(63) : error C2065: 'output' : undeclared identifier

diff_eq3.c(63) : warning C4022: 'memcpy' : pointer mismatch for actual parameter 1

diff_eq3.c(67) : error C2065: 'delete' : undeclared identifier

diff_eq3.c(67) : error C2059: syntax error : ']'

diff_eq3.c(68) : error C2065: 'delete' : undeclared identifier

diff_eq3.c(68) : error C2059: syntax error : ']'

Are you sure you are compiling in C++ and not C (your compiler probably has the option to let you do one or the other)?

ReplyDelete