Setting up a system of equations for ODE45?
4 views (last 30 days)
Show older comments
I had asked this question yesterday but quickly realized i made some invalid assumptions. I need to set up a system of equations to solve with ode45 for the following equations. My goal is to maximize E.
dN/dt = k6*N
dE/dt = (k2*N*S)/(k5+S) - k3*N^2 - k4*E*t
dS/dt = -k7*9
The k values are as follows:
k2 = .21
k3 = 4*10^-8
k4 = .1
k5 = .05
k6 = 1
k7 = 2*10^-7
Any help would be greatly appreciated. I know how to set up for the first and last equations but that middle one messes me up.
1 Comment
Accepted Answer
Walter Roberson
on 6 Dec 2011
The boundary conditions N(0) and S(0) are needed to address this problem. Let N(0) = P and S(0) = Q; then
E(t) = int(-((1/25000000)* P * (1/20-(9/5000000)*_z1+Q) * exp(_z1) + (17/44973545)*_z1 - (21/100)*Q) * P * exp(-(1/20)*t^2 + _z1 + (1/20)*_z1^2) / (1/20-(9/5000000)*_z1+Q), _z1 = 0 .. t)
There does not appear to be an analytical solution for the integral.
_z1 here is a dummy variable for integration purposes.
If we assume that P and Q are real-valued, and that t is non-negative, then by inspection of the numerator, we can see that the numerator is real-valued.
Now examine the denominator, (1/20-(9/5000000)*_z1+Q) where _z1 is the dummy variable over 0 to t.
If Q > -1/20 then this starts out negative and crosses to positive, which implies a division by 0, leading to an undefined integral. The integral is defined, though, if Q < -1/20.
If we arbitrarily set Q = -1 then it becomes possible to probe the conditions under which we can maximize E(t) . Cutting short a bunch of investigation, E(t) goes negative for sufficiently large t. How large is "sufficiently large"? Best case: after t=1.1, when P = 2.763157895*10^6 . That is, for larger or smaller P, E(t) would go negative sooner. Using a different Q < -1/20 would not change the basic behavior, but can change the P that maximizes the expression, and so shift the location of the peak.
I would say that since E(t) will go negative sooner or later (usually sooner!), you cannot really speak of maximizing E(t)
3 Comments
Walter Roberson
on 6 Dec 2011
(Reminder for this discussion: P is the boundary condition N(0) )
I explored this further, and I thought I had encountered a situation where the peak was higher than I expressed above. That turned out not to be the case, though: the P = 2.763157895*10^6 does give the maximum E(t) value, with t about equal 0.8.
What I did find was that as P increases up to the above value, the peak squeezes to lower and lower t values, but the maximum E(t) increases. Conversely, as P approaches 0 from above, the location of the peak E(t) increases; for example at P=1/1000 the peak is slightly to the right of t=21. At P=1E-9, the peak is about 35.75. As P gets smaller, the curve continues to flatten out; at the limit P=0, the curve is flat.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!