How to introduce a maximum value constraint when solving for variables using ode45?

Hello, I'm trying to solve some ode functions as I'm developing a constitutive model for complex fluid and ran into this issue.

My first variable or the "y" in ode45 is a spacially resolved variable "A" along 1D axis, which changes over time. I have the formula for the time derivative for A, so I can solve for the time evolution of, for example, A1 to A10 (the y(1) to y(10) in ode45) representing A at different positions. That's a pretty standard case for using ode45.

Now I want to introduce a second variable which also changes over time at different positions and lets just call it "B" (or B1 to B10 in reality). I have the formula for the time derivative for B, but B also has a critical value constraint B_max. Therefore, the increase of B should stop once it reaches the maximum value B_max(A), which itself is not a constant but has linear relationship with A at the given position. I'm expecting a scenario that A decreases over time, B initially increases but stops at the B_max and later decreases since B_max decreases with A. B can also decrease after it reaches the maximum value when the time derivative becomes negative.

I've been thinking about how to introduce this maximum constraint but couldn't figure it out. The first way to do it is to compare the B value with B_max and if the increment from the previous time step made B > B_max, then we set B = B_max and no longer use the normal formula for dB/dt and just set dB/dt = 0. However, the problem is that while I can check if B > B_max in the odefun dy = f(y, t) passed to ode45, I could not directly change B value in odefun since because B is not passed by reference as the input and the solver automatically stores and updates B value.

Then I thought about modifying the time derivative that if the increment dB+B > B_max, then set dB = B_max - B so the new B won't go beyond B_max. However, the problem is dB = (dB/dt)*dt and the time step dt is not a fixed value and is optimized by the solver, so I don't know what value of dt the solver is going to choose for the next step and therefore can't determine dB/dt.

Finally, I tried setting B as a global variable and updating it myself in the odefun dy = f(y, t) and let the solver just solve for A. I used a global variable to record the t passed into the odefun and calculated the timestep used by the solver in the previous step by comparing the previous time and current time. I can use that time step to update B with the formula for time derivative and also modify B as I want when it reaches the maximum value. With this method, I ran the matlab code and solved the problem without getting error messages. However, I'm getting fluctuating time evolution curves for B which I believe is due the inappropriate choice of time step by the ode45 solver. Since now it is only solving for A and doesn't know about B, the time step is only optimized for A and could be too large for B.

I'm wondering if there is any way in Matlab to solve my problem. It seems to me there is no way around unless I can write a custom ode45 solver myself, which unfortunately is beyond my capability. Any thought is welcomed. At least if I learn it is something that just could not be done in Matlab, then I can think about alternative plans for my project.

Answer 1

[Steven Lord]

Use an events function.