Python code for second order non linear differetial equation

Question

To find the deformations on the quarter of a ring I found the following formulations:

ϵ=(θ′−1/R)∗t/2, θ(s)′′=−F4EI∗cosθ, θ(0)=π/2, θ(L/4)=0

Being L the medimum circumference, R the medium radius, t the thickness, E young modulus, I inertia and s the spatial coordinate following the arc.

I tried writing a code to solve this but the deformations I obtain are only negative while there should be a change in sign at some angle. Below is the code I tried:

``
import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# Costanti note
F = 10  # Forza massima applicata in Newton
E = 70000  # Modulo di Young in MPascal
R_e = 28  # Raggio esterno in millimetri
R_i = 27.5  # Raggio interno in millimetri
t = R_e - R_i  # Spessore dell'anello in millimetri
h = 5  # Altezza della sezione in millimetri
I = 1/12 * h * t**3  # Momento d'inerzia della sezione in mm^4
R = (R_e + R_i) / 2  # Raggio medio in millimetri
L = 2 * np.pi * R  # Circonferenza media in millimetri

# Equazioni differenziali
def equations(s, y):
    theta = np.pi / 2 - s / R  # Calcolo dell'angolo theta per ogni punto s
    theta_prime = y[1]
    dtheta_prime_ds = - (F / (4 * E * I)) * np.cos(theta)
    return np.vstack((theta_prime, dtheta_prime_ds))

# Condizioni al contorno
def bc(ya, yb):
    return np.array([ya[0] - np.pi / 2, yb[0]])

# Intervallo di integrazione
s = np.linspace(0, L / 4, 100)

# Condizioni iniziali stimate
y_init = np.zeros((2, s.size))
y_init[0] = np.pi / 2

# Risoluzione del problema ai valori al contorno
solution = solve_bvp(equations, bc, s, y_init)

# Verifica della soluzione
if not solution.success:
    raise RuntimeError("La soluzione dell'equazione differenziale non ha avuto successo.")

# Calcolo della soluzione per questi valori di s
s_values = np.linspace(0, L / 4, 100)
theta_values = solution.sol(s_values)[0]
theta_prime_values = solution.sol(s_values)[1]

# Stampa dei valori di theta e theta'
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)

# Calcolo della deformazione epsilon
epsilon_values = (theta_prime_values - (1 / R)) * t / 2

# Stampa dei risultati
print("s:", s_values)
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)
print("epsilon(s):", epsilon_values)

# Visualizzazione dei risultati
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(s_values, theta_values)
plt.xlabel('s (mm)')
plt.ylabel('theta(s)')
plt.title('Soluzione theta(s)')

plt.subplot(3, 1, 2)
plt.plot(s_values, theta_prime_values)
plt.xlabel('s (mm)')
plt.ylabel('theta\'(s)')
plt.title('Derivata theta\'(s)')

plt.subplot(3, 1, 3)
plt.plot(s_values, epsilon_values)
plt.xlabel('s (mm)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon(s)')

plt.tight_layout()
plt.show()

# Grafico separato di epsilon vs theta (in gradi)
theta_degrees = np.degrees(theta_values)
plt.figure(figsize=(8, 6))
plt.plot(theta_degrees, epsilon_values)
plt.xlabel('theta (degrees)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon vs theta (degrees)')
plt.grid(True)
plt.show()
``

Answer 1

I don't think there is anything wrong with your code. However, I think the boundary conditions that were specified in the original source paper are the wrong way round. (The paper is a conference paper, not a peer-reviewed journal one.)

Consider the case F=0. In this case θ(s)′′=0, so θ(s)′=constant. That constant can be worked out as the gradient between the boundary conditions. If you took those boundary conditions as θ(0)=π/2, θ(L/4)=0 as you (following the original paper) have done, then θ(s)′=-2π/L=-1/R. The expression for ϵ is then inevitably negative throughout. But worse, you have ZERO FORCE BUT A NON-ZERO STRAIN! That really would be magic!

I suspect one of two things. Scanning through that paper suggests that θ might be being used for two completely different things (or, at least, with different orientations in different areas). More likely, however, is that the boundary conditions are simply the wrong way round. If you try θ(0)=0, θ(L/4)=π/2 (as looks more natural given that θ is supposed to parameterise the length of arc) then ϵ will change sign as you desire. (An added bonus is that ϵ will also be zero when the force F is zero, which is all to the good.)