How to solve ODEs that are a function of the derivatives of the state variables

I don't know how to properly ask this question. Let me try to explain it.

Typical coupled ODEs:

y'(1) = function(state variables)
y'(2) = function(state variables)
y'(3) = etc ...

Can I solve ODEs in MATLAB that look like this?:

y'(1) = function(state variable,y'(x))
y'(2) = function(state variables,y'(x))
y'(3) = etc ...

basically I am asking if you can have y primes' on the right side of the equal sign.

(Edit 3/1/2011) example:

• rho' = f(rho,A,u,u')
• u' = f(rho,u,P')
• h' = f(u, u')
• P' = f(alpha,T',T,rho',rho)
• T' = f(alpha,h')
• alpha'= f(T,rho',rho,alpha)

It is not possible to write the primes' in terms of state variables only. Can MATLAB solve these type of equations? Is there a similar function to Mathematica's NDSolve[]?

[Matt Tearle]

The notation is a bit confusing. Can you give a small example. I think the answer is "yes", but an example would help show how (as well as understand the question correctly)

[Jesse]

Sure Matt (I agree it's confusing)

I want to solve flow through a rocket nozzle using state variables: density (rho), velocity (u), pressure (P), temperature(T), enthalpy (h), and degree of dissociation (alpha). The independent variable is area (basically a 1-D analysis).

Here are the equations

drho/dA = -(rho/A)-(rho/u)du/dA

du/dA = -(1/(rho*u)) dP/dA

dh/dA = -u*du/dA

dp/dA = R( (1+alpha)*(T*drho/dA + rho*dT/dA))

dT/dA = function(dh/dA,dalpha/dA)

dalpha/dA = function(dT/dA,drho/dA)

Every function contains a d/dA of other state variables.

Typically, all the d/dA's are only on the left (I think)

Thank you for your help

Answer 1

[Matt Tearle]

Ah, that actually wasn't what I was thinking (I don't know why, b/c in retrospect your question makes perfect sense). Anyway, that actually makes the answer fairly easy: yes. See ode15i. Rewrite the equations as a vector system F(A,x,x'), where x = [rho; u; h; P; T; alpha], and away you go.

[Jesse]

Matt,

I think you're right. I say think because it looks like the right ode solver, but I can't get it to work. I suspect it failed for two reasons;

1 - I don't know the syntax

2 - I don't know the initial d/dA conditions.

Here is my function

==========================================================

function Z = implicit(A, y,yP)

global RO2 thetaD rhoD

Z=zeros(5,1);

rho = y(1);

u = y(2);

h = y(3);

P = y(4);

T = y(5);

a = .92;

rhoP = yP(1);

uP = yP(2);

hP = yP(3);

PP = yP(4);

TP = yP(5);

z=[rho + ( (rho/A) + (rho/u)*uP); u + (1/(rho*u))*PP; hP + u*uP; PP - RO2*( (1+a)*(T*rhoP + rho*TP) ); TP - ( - (1/RO2)*hP )/(4+a)];

============================================================

I made guess for the initial conditions, but have no idea if they satisfy the equations (0's didn't work)

I get this error:

Warning: Failure at t=1.000000e+00. Unable to meet integration tolerances...

Any feedback would be appreciated!

[Matt Tearle]

Works for me, although I saw a couple of typos in your function (not sure if you copy/pasted, or they're really there).

Last line is z, not Z (so it always returns all zeros!). Also the first two elements of Z start with rho and u instead of rhoP and uP.

Having fixed those, I ran it no problem with:

A0 = 1;

y0 = rand(5,1);

yp0 = fsolve(@(yp) implicit(1,y0,yp),rand(5,1))

[Asolve,Zsolve] = ode15i(@implicit,[1,20],y0,yp0);

plot(Asolve,Zsolve)

Note the use of fsolve to get an initial guess for the derivatives. A and the state variables are fixed, so finding the derivatives amounts to finding yP such that implicit(A,y,yP) = 0.

While we're here, using globals is icky. Pass parameters using anonymous function handles. Make implicit a function of A, y, yP, *and* any extra parameters you need. Then do:

foo = %set value

f = @(A,y,yP) implicit(A,y,yP,foo);

[Asolve,Zsolve] = ode15i(f,[1,20],y0,yp0);

This way, f is a handle to implicit with the value of foo "embedded".

[Jesse]

I have not tried your suggestions yet. Let me first say thank you for your help. Wow! so much detail in your explanations.

I should have know better about globals. I read about people getting reprimanded for this all the time :)

And yes, those were errors! I copied and pasted my code.

[Jesse]

I've said it before, but thank you Matt!

Do you have a link to a page that describes function handles? They don't make much sense to me :-(

[Matt Tearle]

Here are the online doc pages:

Function handles in general

http://www.mathworks.com/help/matlab/ref/function_handle.html

Anonymous function handles

http://www.mathworks.com/help/matlab/matlab_prog/f4-70115.html

And, most useful in this context, using anonymous handles to provide parameters to named functions, for use in things like the ode routines

http://www.mathworks.com/help/matlab/matlab_prog/f4-70115.html#f4-70238

Answer 2

[Jesse]

By the way, final code I used: (for other ode15i users)

===================================

function Z = implicit(A, y,yP)
global RO2 thetaD rhoD

Z=zeros(5,1);

rho = y(1);
u = y(2);
h = y(3);
P = y(4);
T = y(5);
a = .92;

rhoP = yP(1); % rhoP = -((rho/A) + (rho/u)*uP)
uP = yP(2);   %  uP  = -(1/(rho*u))*PP
hP = yP(3);   %  hP  = -(hP + u*uP)
PP = yP(4);   %  PP  = (RO2*(1+a)*(T*rhoP + rho*TP))
TP = yP(5);   %  TP  = (1/RO2)*hP/(4+a)

Z=[ rhoP + ((rho/A) + (rho/u)*uP)       ;
     uP  + (1/(rho*u))*PP               ;
     hP  + (hP + u*uP)                  ;
     PP  - (RO2*(1+a)*(T*rhoP + rho*TP));
     TP  - (1/RO2)*hP/(4+a)             ]

[Jesse]

Please excuse the "globals"

That was naughty of me to use them!

[Matt Tearle]

Thanks for posting that, for the benefit of others. I'd call that an example of "doing the right thing"

(http://www.mathworks.com/company/aboutus/mission_values/)