I encountered a failure on solving deflection curve using scipy.solve_ivp when I used a big length for the beam(pipe). The following is my script.
# -*- coding: utf-8 -*-
#
# Calculate deflection curve under gravity
#
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# round pipe's sectional area
def area_round(D,t):
d = D - 2 * t
return np.pi*(D*D-d*d)/4.0
# area moment of inertia section properties of round pipe
def moi_round(D,t):
d = D - 2 * t
return np.pi/64.0*(D**4-d**4)
# gravity force per unit length for round pipe
def G_round(dens,D,t):
g = 9.8 # gravity coef, N/kg
A = area_round(D,t)
return A*dens*g
# Calculate 2D plane(xy) bending of a beam undering self-gravity
# minus y is the gravity direction
# - Young's modules E
# - area moment of inertia: Iz
# x: coordinate position along x-axis
# u=[y,dydx]: deflection at y direction , first order derivative on x
def deflection(x,u,E,Iz,G,LG):
y,dy_dx = u
MG = np.select([x<=LG,x>LG],[G*(0.5*LG**2+0.5*x**2-LG*x),0]) # distributed force
ddy_ddx = np.power(1.0+dy_dx*dy_dx,3.0/2.0)*MG/E/Iz
return [dy_dx,ddy_ddx]
# solve initial value problem using numerical method
# Given x list, calcuate y and y'
def solve_deflection(E,Iz,Gy,Lg,x_list):
# initial value
y0 = 0.0
dy_dx_x0 = 0.0
sol = solve_ivp(deflection,t_span=[0.0,Lg],y0=[y0,dy_dx_x0],method='RK45',t_eval=x_list,vectorized=True,args=(E,Iz,Gy,Lg),rtol=1.0e-3,atol=1.0e-4)
print("sol:",sol)
return sol.y[0,:],sol.y[1,:]
# test of solving deflection equation
# unit sytem: length mm force:N pressure:N/mm^2/MPa
def test():
dens = 7880.0e-9 # Kg/mm^3
E = 210e3 # elastic modulus, N/mm^2
# round pipe
Dp = 60.0 # mm, outer diameter
tp = 3.0 # mm ,thickness
print("Dp:",Dp,' tp:',tp)
Iz= moi_round(Dp,tp) # section const
Gy = -G_round(dens,Dp,tp) # gravity force per unit
print(" Iz:",Iz,' Gy:',Gy)
# Lg = 18799.1 # mm, this length is OK
Lg = 18799.2 # mm, this length leads to failure
print("Lg=",Lg)
xlist = np.linspace(0.0,Lg,num=100,endpoint=True)
res = solve_deflection(E, Iz, Gy=Gy, Lg=Lg, x_list=xlist)
plt.plot(xlist,res[0])
plt.show()
if __name__ == '__main__':
test()
pass
When the length is 10000mm, it can give correct result as the following graph
But when I changed it to 20000mm, it failed and output the following messages,
sol: message: 'Required step size is less than spacing between numbers.'
nfev: 944
njev: 0
nlu: 0
sol: None
status: -1
success: False
t: array([ 0. , 202.02020202, 404.04040404, 606.06060606,
808.08080808, 1010.1010101 , 1212.12121212, 1414.14141414,
1616.16161616, 1818.18181818, 2020.2020202 , 2222.22222222,
2424.24242424, 2626.26262626, 2828.28282828, 3030.3030303 ,
3232.32323232, 3434.34343434, 3636.36363636, 3838.38383838,
4040.4040404 , 4242.42424242, 4444.44444444, 4646.46464646,
4848.48484848, 5050.50505051, 5252.52525253, 5454.54545455,
5656.56565657, 5858.58585859, 6060.60606061, 6262.62626263,
6464.64646465, 6666.66666667, 6868.68686869, 7070.70707071,
7272.72727273, 7474.74747475, 7676.76767677, 7878.78787879,
8080.80808081, 8282.82828283, 8484.84848485, 8686.86868687,
8888.88888889])
t_events: None
y: array([[ 0.00000000e+00, -3.66155907e+00, -1.45618988e+01,
-3.25951385e+01, -5.76794150e+01, -8.97565181e+01,
-1.28794030e+02, -1.74786625e+02, -2.27730361e+02,
-2.87640987e+02, -3.54556806e+02, -4.28538678e+02,
-5.09670018e+02, -5.98056798e+02, -6.93827544e+02,
-7.97133338e+02, -9.08147821e+02, -1.02706718e+03,
-1.15411018e+03, -1.28951811e+03, -1.43356488e+03,
-1.58659673e+03, -1.74895654e+03, -1.92104924e+03,
-2.10335332e+03, -2.29642078e+03, -2.50087718e+03,
-2.71742163e+03, -2.94684282e+03, -3.19021861e+03,
-3.44843893e+03, -3.72265946e+03, -4.01447603e+03,
-4.32592470e+03, -4.65948174e+03, -5.01806360e+03,
-5.40502733e+03, -5.82494500e+03, -6.28452879e+03,
-6.79299154e+03, -7.36408180e+03, -8.02047137e+03,
-8.80421033e+03, -9.81478548e+03, -1.15113253e+04],
[ 0.00000000e+00, -3.61399978e-02, -7.16869806e-02,
-1.06770973e-01, -1.41512819e-01, -1.76024795e-01,
-2.10415154e-01, -2.44799321e-01, -2.79258778e-01,
-3.13875788e-01, -3.48738932e-01, -3.83943102e-01,
-4.19589504e-01, -4.55785660e-01, -4.92645404e-01,
-5.30288885e-01, -5.68842566e-01, -6.08439223e-01,
-6.49217947e-01, -6.91324142e-01, -7.34921284e-01,
-7.80235186e-01, -8.27461362e-01, -8.76842769e-01,
-9.28683414e-01, -9.83348358e-01, -1.04126372e+00,
-1.10291665e+00, -1.16888799e+00, -1.24022057e+00,
-1.31734604e+00, -1.40101902e+00, -1.49267854e+00,
-1.59444811e+00, -1.70913565e+00, -1.84023352e+00,
-1.99191927e+00, -2.17052401e+00, -2.38680965e+00,
-2.65797898e+00, -3.01398821e+00, -3.51651158e+00,
-4.31409773e+00, -5.92506745e+00, -1.41760305e+01]])
y_events: None
Traceback (most recent call last):
File "mwe.py", line 82, in <module>
test()
File "mwe.py", line 76, in test
plt.plot(xlist,res[0])
File "\lib\site-packages\matplotlib\pyplot.py", line 3019, in plot
return gca().plot(
File "\lib\site-packages\matplotlib\axes\_axes.py", line 1605, in plot
lines = [*self._get_lines(*args, data=data, **kwargs)]
File "\lib\site-packages\matplotlib\axes\_base.py", line 315, in __call__
yield from self._plot_args(this, kwargs)
File "\lib\site-packages\matplotlib\axes\_base.py", line 501, in _plot_args
raise ValueError(f"x and y must have same first dimension, but "
ValueError: x and y must have same first dimension, but have shapes (100,) and (45,)
I have read similar question on solve_ivp error: 'Required step size is less than spacing between numbers.' But I can not understand it.
I also tried different atol value and can not solve the issue. I have struggled on this for several days.
Anyone can help me out? Thank you.
