Solving multiple coupled ODEs

Python Help
1

Hi, I’m fairly new to more advanced python code but I’m trying to solve three coupled ODEs using odeint. I keep getting a message and then an output of wacky looking graphs:

Message
RuntimeWarning: overflow encountered in scalar multiply
dmdt = (0.413 * A * gamma * (rho * math.exp(-h/H)) * v**3)/HVAP
ODEintWarning: Illegal input detected (internal error). Run with full_output = 1 to get quantitative information.
x = odeint(odes,x0,t)
RuntimeWarning: overflow encountered in power
return np.power(self.base, values)

code

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import math

cd = 2
rho = 1.12e-08 # kg/m^3
A = math.pi * (5.811812)**2   # m^2
gamma = 1
HVAP = 80800000              # J/kg

theta = 45 *math.pi/180

H = 85000 # scale height (m)

# rho_atm = (rho * math.exp(-h/H))

# v_0 = 18.14297*10**5             
# m_0 = 3025741
# h_0 = 140000   

def odes(x,t):
    v = x[0]
    m = x[1]
    h = x[2]

    dvdt = (-1/2 * cd * (rho * math.exp(-h/H)) * v**2 * A)/m
    dmdt = (0.413 * A * gamma * (rho * math.exp(-h/H)) * v**3)/HVAP
    dhdt = -v*math.sin(theta)

    return [dvdt, dmdt, dhdt]

x0 = [18.14297*10**5, 3025741, 140000]

# print(odes(x0,0))

t = np.linspace(0,15,1000)
x = odeint(odes,x0,t)

v = x[:,0]
m = x[:,1]
h = x[:,2]

plt.plot(t,v)
plt.show()

plt.semilogy(t,m)
plt.show()

plt.plot(t,h)
plt.show()
2

Could you, please, include all code to reproduce error? Provided part include a lot of parameters, it’s not clear also which initial values you are using.

4

The ODE has a singularity:

    v = x[0]
    m = x[1]
dvdt = ( ...  )/m

More analysis is needed (e.g. asymptotics) to decide what to do in the (v,m,h) phase space near m=0, and to ensure the domain being numerically solved on is chosen s.t. ideally the IVP’s well-posed, or at least that the solver will converge to the desired solution.

Also a note of caution. Ploughing onwards determinedly, and solving equations for their own sake is all very well (afterall, doing so is fun and often results in pretty graphs). But if the solution is going to be relied upon or treated as meaningful, singularities ought to prompt setting aside a moment or two for a few extra checks, and even confirming the interpretation of what the singularity represents, according to the principles and assumptions being used the model, is as intended, and is being treated by the model in the most appropriate way, and does not break any laws of physics.

1 Like
5

Generally when you post code to a forum like this it is better to simplify the code a bit especially to make it self-contained. The example you show could be small enough that I could see that the code is safe to run and copy paste it to try it out. However the fact that I would first need to install openpyxl, make a spreadsheet, etc means that I am not going to do that. It would be better to put the numbers from the spreadsheet directly in the code.

3 Likes
6

My guess as to the problem is that your equations are unstable and h becomes negative and then exp(-h/H) becomes a huge number. I don’t know whether that is because the physical system is actually unstable or if it is because the code or parameters are incorrect or it is down to some numerical inaccuracy though.

It is possible that someone could give a better answer if you show self-contained code that others can try.

1 Like
7

Hello,

as with any problem, you should always try and break it down into its smallest elements.

Let’s start here - the traceback error. Note that it references the following equation:

RuntimeWarning: overflow encountered in scalar multiply

dmdt = (0.413 * A * gamma * (rho * math.exp(-h/H)) * v**3)/HVAP

You should write a test script where you test this equation in isolation. If you need help,
instead of providing this:

A = math.pi * (sheet[f"T{row}"].value)**2

Provide A with the actual value and NOT expect us to read it from some worksheet. In other words, you should provide the value up front. Example:

A = math.pi * (value_here ** 2)  # provide the actual value

What you want to verify is that the equation works with the values that you are providing.

cd = 2
rho = 1.12e-08
gamma = 1
H = 85000

# Provide the values as opposed to having to read them from some spreadsheet
# Provide the raw values
A = math.pi * (sheet[f"T{row}"].value)**2
HVAP = sheet[f"S{row}"].value

dmdt = (0.413 * A * gamma * (rho * math.exp(-h/H)) * v**3)/HVAP

If we use a backpacking analogy, prior to going on said excursion, you first want to make sure that you have all of the items on your check off list: bear spray, water filter, fire starter, …, etc. You don’t wait and verify this until you’re up on some hill somewhere. You do this beforehand. The same applies here. As you’re building your script, you want to make sure and verify that the equations that you’re creating are being tested with the expected hypothetical values and that they’re working as expected prior to adding additional code. This way, if the equation does not work, you’ll know exactly where the issue is coming from. If you wait until you develop an application to test your script, it will make it more difficult to isolate the source of the issue. You should be continously testing as you’re adding additional lines of script.

Do this for the other two equations.

8

The way that you have your function defined, it returns three different equation results in a list. Is this correct / allowable? Are you supposed to return one only?

Reference this link:

If you want similar results, you may have to define three separate functions. One for each of the three equations.

Begin by working on one equation first (verifying and plotting it) so that you know the script works in its simplist form. Once having done so, then add the other two equations / functions.

9

Hi everyone,

Thank you all for the feedback, apologies for the code not being self-contained code. I understand how I could’ve made it easier for you guys to test it yourself, I’ve never posted or used a python forum so I’m not really familiar with the etiquette, but I’ll definitely take your feedback into account next time!

As for the code, I will try testing the equations separately for any singularities and/or other problems stated in your replies.

Edit: I have updated the original post to have the code with all the values, not longer having to use the excel spreedsheet :slight_smile:

10

All this in no way special just for python forums. You expected that people will invest time to help you for free.

Did you check that problem can be reproduced with that code on your system?

I tried to run you example (slightly changed, just to get only a single plot, see below). I didn’t check if this solution is not just a numeric garbage, but at least I see no errors:

xxx
xxx640×480 22.1 KB

Maybe error happens for some specific numpy/scipy versions?

>>> numpy.__version__
'2.2.3'
>>> scipy.__version__
'1.15.2'
>>> 

# xxx.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import math

cd = 2
rho = 1.12e-08 # kg/m^3
A = math.pi * (5.811812)**2   # m^2
gamma = 1
HVAP = 80800000              # J/kg

theta = 45 *math.pi/180

H = 85000 # scale height (m)

def odes(x,t):
    v = x[0]
    m = x[1]
    h = x[2]

    dvdt = (-1/2 * cd * (rho * math.exp(-h/H)) * v**2 * A)/m
    dmdt = (0.413 * A * gamma * (rho * math.exp(-h/H)) * v**3)/HVAP
    dhdt = -v*math.sin(theta)

    return [dvdt, dmdt, dhdt]

x0 = [18.14297*10**5, 3025741, 140000]

t = np.linspace(0,15,1000)
x = odeint(odes,x0,t)

v = x[:,0]
m = x[:,1]
h = x[:,2]

ax1 = plt.subplot(212)
ax1.plot(t,v)
ax2 = plt.subplot(221)
ax2.semilogy(t,m)
ax3 = plt.subplot(222)
ax3.plot(t,h)
plt.show()
11

To make your graphs a bit more readable, you can modify the plotting script like this. This way you know which graph is which and what is being graphed with respect to time.

plt.plot(t,v)
plt.xlabel('time (s)', fontweight='bold')
plt.ylabel('dv/dt', fontweight='bold')
plt.title('dv/dt vs. time', fontweight='bold')
plt.show()

plt.semilogy(t,m)
plt.xlabel('time (s)', fontweight='bold')
plt.ylabel('dm/dt', fontweight='bold')
plt.title('dm/dt vs. time', fontweight='bold')
plt.show()

plt.plot(t,h)
plt.xlabel('time (s)', fontweight='bold')
plt.ylabel('dh/dt', fontweight='bold')
plt.title('dh/dt vs. time', fontweight='bold')
plt.show()

I noticed that the graph dm/dt vs. time has extremely large y axis values. The maximum value is: 1 x 10^97!. This is generally not practical. In code, aside from there being potentially syntax errors, there may also be mathematical errors that will not show up as traceback errors. This is something that, from reasoning, or experience, you will have to determine if the outcome makes sense. Coincidentally, I was watching this video yesterday.

Here is how to obtain the value of 52!:

import math

value = math.factorial(52)
sci_notation = "{:e}".format(value)
print(sci_notation)

>>> 8.065818e+67  # Value of 52! in scientific notation

# This is much smaller (by trillions trilions ...) than 1 x 10^97

This is why I am hesitant to believe that this graph in particular is practical. Check your math and/or the values that you’re passing in.

Q: Why is it now working fine and not before:
Next, since by plugging in the raw values you are no longer getting syntax errors, the next thing that you can do is verify that you are reading the values from the spreadsheet correctly. You can read the values from the spreadsheet and print them out to verify that they match the values from the spreadsheet. This should be done independently from this script.

Once you have verified the values (after fixing your code) are correct, then proceed as per your original script.

Update:

# x0 = [18.14297*10**5, 3025741, 140000]
x0 = [18.14297*10**2, 3025741, 140000]  # I replaced **5 with **2 in x[0]

This fixes the out of bound y-axis values. Does this supposed to be to the power ot 2 vs to the power of 5? Maybe this is not the root cause but you have to look into why the results are unbounded.

12

It turns out I went about this the hard way :confused: (apologies if I cant explain this very well)

The data that I was using from the excel spreadsheet was the output of a Monte Carlo program I made. It turns out the functions I needed plotted (mass, velocity and height) were updated over a for loop, which acted as the coupled system of equations as each value was updated and affected the next iteration.

Long story short, I took these values from each iteration into a list and plotted them.
I would say it was a waste of time, but I enjoyed learning how to integrate coupled ODEs.
Thank you to those who replied to help.

height
height640×480 15.9 KB

MassLoss
MassLoss640×480 16.1 KB

Velocity
Velocity640×480 15.8 KB