Setting up a system of equations for ODE45?

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.

[Walter Roberson]

Is E(0)=0 as in the previous iteration of the question?

Answer 1

[Walter Roberson]

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)

[Donald]

Actually this is perfect thank you. I expected it to have a maximum around t=1 or so, so this works out like i had expected.

[Walter Roberson]

(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.

[Donald]

Wow. I can't thank you enough for this. You have helped me so much. Thank you.