Implementing a first order forward difference scheme in MATLAB

Hi everyone,

I am trying to solve the first order differential equation on the figure below:

So that I can reproduce the curves below:

As you can see from the curves the specific heat has only two values 1.66 and 1.4. In the X equation, all the variables except x are constant. x can be taken to increase from 0 to 10. I want to produce the variation of M with X. Here they used first order forward difference, I'd like to know how I could do this using MATLAB. I'd appreciate any form of assistance. Regards Charles

[Lukas Bystricky]

MATLAB has their own ODE solvers built in.

Alternatively, if you're interested in replicating that graph exactly (using first order methods), you should look up the forward euler method. It's very simple and requires little more than a single for loop.

Answer 1

[Torsten]

fun=@(X,M)(((2+(gamma-1)*M^2)/2)^(-1/(gamma-1)))/((M-1)*(2+(gamma-1)*M)/(2+(gamma-1)*M^2));
gamma=1.4;
[X1,M1]=ode45(fun,[0 10],1.001);
gamma=1.66;
[X2,M2]=ode45(fun,[0 10],1.001];
plot(X1,M1(:,1),X2,M2(:,1));

Best wishes

Torsten.

[Charles Kubeka]

Thanks for the quick response guys Torsten I tried to run your script and I got the following error: Error: Unbalanced or unexpected parenthesis or bracket. I fixed that as follows:

gamma=1.4; fun=@(X,M)((2+(gamma-1)*M^2/2)^(-1/(gamma-1)))/((M-1)*(2+(gamma-1)*M)/(2+(gamma-1)*M^2)); [X1,M1]=ode45(fun,[0 10],0); gamma=1.66; [X2,M2]=ode45(fun,[0 10],0); plot(X1,M1(:,1),X2,M2(:,1));

But now I get this:

The form of the curve looks different from what they got in the figure. I Appreciate your assistance.

[Torsten]

Should be corrected, but better check again.

Best wishes

Torsten.

[Lukas Bystricky]

The brackets are a little incorrect:

 (2+(gamma-1)*M^2/2)

should be

   ((2+(gamma-1)*M^2)/2)

for instance.

[Torsten]

I changed the code ; the initial condition is M(X=0)=1 instead of M(X=0)=0.

Best wishes

Torsten.

[Charles Kubeka]

I made some corrections in the brackets now I get this:

The curve seems like a mirror of the original. The M values are supposed to be positive. Thanks again.

[Lukas Bystricky]

You didn't change the initial condition. M(0) = 1 (or thereabouts) from the original graph.

[Torsten]

You took the code from above with the change in the initial condition for M ?

Best wishes

Torsten.

[Charles Kubeka]

I just did and now I get this:

This is much much better, thanks a lot guys. Do you have an idea why the other curve is not showing or are they identical? Also how can I extract the values of M for each X value?