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

2011 2 25
2 个回答
回答已采纳
3 次查看（30 天）

0 个投票

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[]?

2 个评论
显示 无 隐藏 无

Matt Tearle
Matt Tearle 2011-2-25
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
Jesse 2011-3-1
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

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Matt Tearle
Matt Tearle 2011-3-2

0 个投票

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.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Jesse
Jesse 2011-3-2
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
Matt Tearle 2011-3-2
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
Jesse 2011-3-2
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
Jesse 2011-3-2
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
Matt Tearle 2011-3-2
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

请先登录，再进行评论。

更多回答（1 个）

Jesse
Jesse 2011-3-2

1 个投票

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) ]

2 个评论
显示 无 隐藏 无

Jesse
Jesse 2011-3-2
Please excuse the "globals"
That was naughty of me to use them!
Matt Tearle
Matt Tearle 2011-3-2
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/)

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by