Combining multiple coupled 1st order equations in python

Question

I'm having serious troubles with solving translating 3 coupled differential equations into python.

The 3 DE's stem from a 4th order DE used to calculate the bending moment of an underwater pipeline that has been substituted/transformed into 1st order equations and are coupled to each other.

4th order:

$\frac{d^4 w(x)}{dx^4} - \frac{Q_{grnd}+ Q_{fric}}{EI} w(x) = 0$

3x 1st order:

$w_1(x_1) = A x_1^3 + B x_1^2 + C x_1 + D + \frac{(q_{subm} x_1^4)}{(24 EI)}$
$w_2(x_2) = E x_2^3 + F x_2^2 + G x_2 + H + \frac{(q_{subm} x_2^4)}{(24 EI)}$ $w_3(x_3) = K e^{(\sqrt2 \beta x_3)}+ L e^{(-\sqrt2 \beta x_3)} + M \cos(\sqrt2 \beta x_3)+N \sin(\sqrt2 \beta x_3)- \frac{(q_{subm})}{(Q_{grnd}+Q_{fric})} $

where $0 \leq x_1 \leq I_1$, $0 \leq x_2 \leq I_2$ and $0 \leq x_3 \leq I_3$

The system has 15 unknown constraints and 15 boundary constraints conditions for which i want to solve.

Boundary conditions

1. $w_1(0) = 0$ (deflection at $x_1$ = 0 is 0), D=0
2. $w_1'(0) = 0$ (slope at $x_1$ = 0 is 0), C=0
3. $w_1''(0) = 0$ (bending moment at $x_1$ = 0 is 0), B=0
4. $w_1(I_1) = -p$
5. $w_2(0) = p$, H=p
6. $w_1'(I_1) = w_2'(0)$
7. $w_1''(I_1) = w_2''(0)$
8. $w_2(I_2) = 0$
9. $w_3(0) = 0$
10. $w_2'(I_2) = w_3'(0)$
11. $w_2''(I_2) = w_3''(0)$
12. $w_2'''(I_2) = w_3'''(0)$
13. $w_3(I_3) = d$
14. $w_3(I_3)' = 0$
15. $w_3(I_3)'' = 0$

I've looked online and tried some codes myself, but until now I unfortunately haven't got any useful results. An example of my code is given below, where 15 start / 'guess' values are given for the unknown variables.

I1 = 46.483
I2 = 5.916
I3 = 21.90



A = -3.858*10**(-5)
B = 0
C = 0
D = 0
E = -1.928*10**(-4)
F = 4.270*10**(-3)
G = 0.049
H = p
K = 0.052
L = 0.1
M = -0.021
N = 0.787


def f(w, t):
     x1 = w[0]
     x2 = w[1]
     x3 = w[2]
     # the model equations (see Munz et al. 2009)

     F0 = A*x1**3 + B*x1**2 + C*x1 + D + (q_subm*x1**4)/(24*EI)
     F1 = E*x2**3 + F*x2**2 + G*x2 + H + (q_subm*x2**4)/(24*EI)
     F2 = K*np.exp(np.sqrt(2)*beta*x3)+L*np.exp(-np.sqrt(2)*beta*x3)+M*np.cos(np.sqrt(2)*beta*x3)+N*np.sin(np.sqrt(2)*beta*x3)-(q_subm/(Q_grnd+Q_fric)) 
     return [F0, F1, F2]

w0 = [0, 0, 0]
t = np.linspace(0, 80, 1000) 
soln = odeint(f,w0,t)

w1 = soln[:,0]
w2 = soln[:,1]
w3 = soln[:,2]

plt.figure()
plt.plot(t, w1, label='w1')
plt.plot(t, w2, label='w2')
plt.plot(t, w3, label='w3')
plt.xlabel('x')
plt.ylabel('deflection')
plt.title('def')
plt.legend(loc=0)

Could anybody please be so kind to assist me?

With kind regards,

Ronald

Answer 1

I think there are some additional assumptions in your equations. A fourth order equation reduces to four first order equations. Given the fourth order equation $$\frac{d^4 w(x)}{dx^4} - \frac{Q_{grnd}+ Q_{fric}}{EI} w(x) = 0$$ the corresponding set of first order equations with $w(x)=\left(\begin{array}{c}w_1(x)\\w_2(x)\\w_3(x)\\w_4(x)\end{array}\right)$ reads: $$\frac{dw}{dx}=\left(\begin{array}{ccc}0&1&0&0\\0&0&1&0\\0&0&0&1\\\frac{Q_{grnd}+ Q_{fric}}{EI}&0&0&0\end{array}\right)w$$ You can solve this numerical and the solution will be in $w_1(x)$.

However, there is no need for a numerical solution, if $S=\frac{Q_{grnd}+ Q_{fric}}{EI}$ is not dependent on $x$. Given the matrix $M$ of the first order system, you simply want to calculate the eigenvalues and eigenvector of $M$ to obtain the solution of the system. The solution is a linear combination of the four eigenvectors $w(x)=\sum_{i=1}^4a_i \phi_i(x)$. The unknown $a_i$ are given by the boundary conditions.

Edit

Thanks for the clarification. The domain $[x_0,x_3]$ with $x_i=\sum_j^i I_j$ of the fourth order ODE is decomposed in three subdomains, $[x_0,x_1],[x_1,x_2],[x_2,x_3]$.

You can do two things, either implement the given ode with $x$-dependent $S$ together with four boundary conditions and obtain the solution on the entire domain.

Or you demand that the solutions in different subdomains are matched along their common boundary by requiring that $w(x)$ and a finite number of its derivatives are equal along the interface. The number of matching conditions is the number that is necessary to give a unique solution in each subdomain, here four.

So possible choices for the interface conditions are e.g. symmetric $x_j\in\{x_1,x_2\}$ and $i\in\{1,2\}$ $$\begin{eqnarray}w_1(x_0)&=&\text{left bc1}\\w'_1(x_0)&=&\text{left bc2}\\w_i(x_j)&=&w_j(x_j)\\ \frac{dw_i}{dx}(x_j)&=&\frac{dw_j}{dx}(x_j)\\w_3(x_3)&=&\text{right bc1}\\w'_3(x_3)&=&\text{right bc2}\end{eqnarray}$$ or asymmetric: $$\begin{eqnarray}w_1(x_0)&=&\text{left bc1}\\w'_1(x_0)&=&\text{left bc2}\\w''_1(x_0)&=&\text{left bc3}\\w_1(x_1)&=&w_2(x_1)\\\\w_2(x_2)&=&w_3(x_2)\\w'_2(x_2)&=&w'_3(x_2)\\w''_2(x_2)&=&w''_3(x_2)\\w_3(x_3)&=&\text{right bc}\end{eqnarray}$$ As you can see, there are always four conditions for each $w_i$. You state six conditions for $w_1$, seven for $w_2$ and seven for $w_3$. These are simply too many and cannot give you a solution.

So you need to think which of your stated boundary conditions are the ones you want, with this equipped you need to implement a system of linear equations for the coefficients $A \dots N$ defined by the conditions, given that the solutions $w_1$ to $w_3$ are valid for each domain.

My guess is, that you want to keep $w_1(0)=0$,$w'_1(0)=0$,$w''_1(0)=0$ and $w_3(x_3)=d$ and use it together with the asymmetric conditions.