0
$\begingroup$

I'm currently trying to solve a problem using numerical methods. The set-up is rather long, so I apologize in advance...

TL;DR: My solutions change depending on how big my steps are and I don't know why.


Set-up

The first part of this problem can be described by a Lagrangian $L(s,\theta)$, which is given by $$L = \frac{1}{24}(L^2\dot\theta^2+4b^2\dot\theta^2+24s\dot\theta^2+12\dot s^2+ 12 b\dot s \dot\theta+24gs\sin\theta -12bg\cos\theta).$$ This Lagragian gives, by means of the Euler-Lagrange equations, $$ \begin{align*} \dot\theta^2 +g \sin\theta &=\ddot s + \frac{1}{2}b\ddot\theta\\ \left(\frac{1}{12}L^2+\frac{1}{3}b^2\right)\ddot\theta+\frac{1}{2}b\ddot s &= \frac{1}{2}bg\sin\theta-2\dot s\dot\theta-gs\cos\theta, \end{align*} $$ where you can assume $L=0.1,\,b=0.02$ and $g=9.81$. I also have the initial conditions given by: $$s(0)=0.1,\quad \dot s(0)= 0.3,\quad u(0)= 0,\quad \dot u(0)= 0.$$ The generalized coordinates $s,\theta$ and the real trajectory are connected over $$x = \begin{pmatrix}s\cos\theta + \frac{b}{2}\sin\theta\\ -s\sin\theta +\frac{b}{2}\cos\theta\end{pmatrix}.$$

The first part of the problem is to solve this set of differential equation up to the time $t^*$ when the condition $\ddot x(t^*) = 0$ is fulfilled, or in terms of $s$ and $\theta$ until we have $$-2 \sin\theta(t^*) \dot s(t^*)\dot\theta(t^*) - \cos\theta(t^*) s(t^*)\dot\theta^2(t^*) - \frac{1}{2}b\sin\theta(t^*) \dot\theta^2(t^*) + \cos\theta(t^*) \ddot s(t^*) + \frac{1}{2} b\ddot\theta \cos\theta(t^*) -s(t^*) \ddot\theta(t^*) \sin\theta(t^*) \overset{!}{=} 0.$$ After this time $t^*$ the new trajectory is given by

$$\vec{x}= \begin{pmatrix}x\\ - \left[x\tan(|\alpha|)+\frac{gx^2}{2v_0\cos^2(|\alpha|)}\right]\end{pmatrix} \quad\text{and}\quad \theta(t)= \dot \theta(t^*) \cdot t + \theta(t^*)$$ where $\alpha = \theta(t^*)$ and $v_0 = \dot s(t^*)$.


My attempt to solve this goes as follows: First defined
$$y = \begin{bmatrix}s\\\dot s\\ \theta\\ \dot\theta\end{bmatrix} \Longrightarrow \dot y = \begin{bmatrix}\dot s\\\ddot s\\\dot \theta\\ \ddot\theta\end{bmatrix}$$ I solved the Euler-Lagrange eq. for $\ddot \theta$ and $\ddot s$ and inserted them into $\dot y$ which gives rhs

def rhs(y, t, g=9.81, L=0.01):
    return np.array([y[1], 
                    (1/(b**2 + L**2 + 24.0**y[0]**2))*
                    (6.0*b*g*np.cos(y[2])*y[0] - 7.0 *b**2 *g*np.sin(y[2]) - g* L*2.0 *np.sin(y[2]) - 24.0 *g *y[0]**2*np.sin(y[2]) 
                        - 24.0*b* y[0]* y[1]* y[3] - 8.0* b**2* y[0] * y[3]**2 - 2.0* L**2* y[0]* y[3]**2 - 48.0*y[0]**3* y[3]**2), 
                    y[3], 
                    -12.0*(-g *np.cos(y[2]) * y[0] + b* g* np.sin(y[2]) + 4* y[0] * y[1]* y[3] + b *y[0] *y[3]**2)/(b**2 + L**2 + 24 *y[0]**2)
                    ])

I then used scipy's odeintalgorithm to solve the problem in the follwoing way

T = np.linspace(0,10,10000)
b = 0.02
lstart = 0.02
ystart = np.array([np.sqrt(lstart**2+ (b/2.0)**2), v, 0.0, 0.0])
y, data = odeint(rhs, ystart, T, full_output = True, args=(g, L))

This worked out and plotting the solutions gives results that I expected.

The next thing to do was to figure out $t^*$. For this I interpolated the solutions and inserted them into the defining condition to find $t^*$. The following function is supposed to return the index of T at which this happens

def find_tstar_index(T, y):
    # Find point where bread leaves table
    s_spl = UnivariateSpline(T, y[:, 0], s=0, k=4)
    u_spl = UnivariateSpline(T, y[:, 2], s=0, k=4)



    s_spl_d = s_spl.derivative(n=1)
    s_spl_2d = s_spl.derivative(n=2)

    u_spl_d = u_spl.derivative(n=1)
    u_spl_2d = u_spl.derivative(n=2)

    searching_for_zeros = -2.0*np.sin(u_spl(T))*s_spl_d(T)*u_spl_d(T) \
                        - np.cos(u_spl(T))* s_spl(T)* u_spl_d(T)**2 \
                        -1/2.0* b* np.sin(u_spl(T))*u_spl_d(T)**2  \
                        + np.cos(u_spl(T))*s_spl_2d(T)\
                        + 1/2.0* b* np.cos(u_spl(T))* u_spl_2d(T) \
                        - s_spl(T)*np.sin(u_spl(T))*u_spl_2d(T)

    tstar_index = np.where(searching_for_zeros < 1e-25)[0][0]
    return tstar_index

tstar_index = find_tstar_index(T, y)

The next step would be to set up $\vec{x}$ and $\theta (t)$ using $t^*$

The actual problem starts here: I now definie the follwoing function

angle = lambda t: y[tstar_index, 3] * t + y[tstar_index, 2]

The thing is now that depending on the T defined in the beginning of my code the function angle changes significantly. For example

T_A = np.linspace(0,1,10000)
T_A = np.linspace(0,10,10000)

give completely different results for the intervall [0,1] when plugged into the respective angle function... There has to be an error in my logic but after hours of trying to figure out what the problem is I'm kind of lost...


Complete Code

I'll post down here my complete code so that you can try it out yourself, but up to line 142 everything seems to work fine (T_fall is the same array regardless of what T is...)

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

g = 9.81
L = 0.1
#T = np.linspace(0,10,10000)
b = 0.02
lstart = 0.02
table_hight = 1.3

vs = [0.7]

def find_tstar_index(T, y):
    # Find point where bread leaves table
    s_spl = UnivariateSpline(T, y[:, 0], s=0, k=4)
    u_spl = UnivariateSpline(T, y[:, 2], s=0, k=4)



    s_spl_d = s_spl.derivative(n=1)
    s_spl_2d = s_spl.derivative(n=2)

    u_spl_d = u_spl.derivative(n=1)
    u_spl_2d = u_spl.derivative(n=2)

    searching_for_zeros = -2.0*np.sin(u_spl(T))*s_spl_d(T)*u_spl_d(T) \
                        - np.cos(u_spl(T))* s_spl(T)* u_spl_d(T)**2 \
                        -1/2.0* b* np.sin(u_spl(T))*u_spl_d(T)**2  \
                        + np.cos(u_spl(T))*s_spl_2d(T)\
                        + 1/2.0* b* np.cos(u_spl(T))* u_spl_2d(T) \
                        - s_spl(T)*np.sin(u_spl(T))*u_spl_2d(T)

    tstar_index = np.where(searching_for_zeros < 1e-25)[0][0]
    return tstar_index

def wurfparabel (x, v, beta):
    b = np.abs(beta)
    return - (x * np.tan(b) + (g * x**2)/(2* v**2* np.cos(b)**2))

def draw_bread(x_p, u, col):

    x_b = x_p + np.array([-b/2.0 * np.sin(u), -b/2.0 * np.cos(u)])
    x_left_down = x_b + np.array([-L/2.0 * np.cos(u), L/2.0 * np.sin(u)])
    x_left_up = x_left_down + np.array([b * np.sin(u), b * np.cos(u)])
    x_right_up = x_left_up + np.array([L * np.cos(u), -L * np.sin(u)])
    x_right_down = x_right_up + np.array([- b * np.sin(u), - b * np.cos(u)])

    plt.plot([x_b[0], x_left_down[0], x_left_up[0], x_right_up[0], x_right_down[0], x_b[0]], 
        [x_b[1], x_left_down[1], x_left_up[1], x_right_up[1], x_right_down[1], x_b[1]], col)


    pass

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(np.abs(array) - value)).argmin()
    return idx

def x_p(y):
    return np.array([y[:, 0] * np.cos(y[:, 2]) + b/2.0 * np.sin(y[:, 2]), 
                    - y[:, 0] * np.sin(y[:, 2] + b/2.0 * np.cos(y[:, 2]))
                    ])

plt.figure()
def test(T):

    def bread_falling_trajectory(v):
        """
        Function that calculates the trajectory of a piece of toast falling from a table with the initial velocity v.
        """

        # Define starting values. The  first value is chosen in such a way that the Center of mass is of the same hight as the zero point.
        ystart = np.array([np.sqrt(lstart**2+ (b/2.0)**2), v, 0.0, 0.0])

        x_sol_start = np.array([ystart[0] * np.cos(ystart[2]) + b/2.0 * np.sin(ystart[2]), 
                        - ystart[0] * np.sin(ystart[2] + b/2.0 * np.cos(ystart[2]))
                        ])

        # Define the right hand side of the differential equations.
        def rhs(y, t, g, L):
            return np.array([y[1], 
                            (1/(b**2 + L**2 + 24.0**y[0]**2))*
                            (6.0*b*g*np.cos(y[2])*y[0] - 7.0 *b**2 *g*np.sin(y[2]) - g* L*2.0 *np.sin(y[2]) - 24.0 *g *y[0]**2*np.sin(y[2]) 
                                - 24.0*b* y[0]* y[1]* y[3] - 8.0* b**2* y[0] * y[3]**2 - 2.0* L**2* y[0]* y[3]**2 - 48.0*y[0]**3* y[3]**2), 
                            y[3], 
                            -12.0*(-g *np.cos(y[2]) * y[0] + b* g* np.sin(y[2]) + 4* y[0] * y[1]* y[3] + b *y[0] *y[3]**2)/(b**2 + L**2 + 24 *y[0]**2)
                            ])

        # Solve the ODE
        y, data = odeint(rhs, ystart, T, full_output = True, args=(g, L))

        # Get the trajectory
        x_sol = x_p(y)

        #Find the index of x_sol at which the bread leaves the table.
        tstar_index = find_tstar_index(T, y)

        # Lösung ab Freier Fall
        # EQM is given by \ddot theta = 0, solving this and inserting the boundry conditions leads to winkel
        winkel = lambda zeit: y[tstar_index, 3] * zeit + y[tstar_index, 2]
        print(y[tstar_index, 3], y[tstar_index, 2] )

        # for the spatial part we essentially have a free falling rod, so we solve the eom and construct trajectory
        v0 = 1/2.0 * np.sqrt(4*y[tstar_index,1]**2 +4*b*y[tstar_index,1]* y[tstar_index,3] + (b**2 + 4*y[tstar_index,0]**2)*y[tstar_index,3]**2)
        beta = y[tstar_index, 2]
        trajectory_ = lambda x: wurfparabel(x, v0, beta)
        y_0 = np.abs(x_p(y)[1][tstar_index]-trajectory_(x_p(y)[0][tstar_index]))

        trajectory = lambda x: trajectory_(x)  + y_0


        # Creat a plot of the results 
        x_free_fall = np.linspace(x_sol[0][tstar_index], 0.45, N)

        sol_dict = {'x_sol_start': x_sol_start, 'table_traje': x_sol,
                    'x_free_fall': x_free_fall, 'free_fall_traj': trajectory(x_free_fall), 'y': y}

        return sol_dict, tstar_index, winkel



    for k, v in enumerate(vs):
        dicct, index, winkel = bread_falling_trajectory(v)
        # Plotting all the results
        #plt.subplot(3,1,k+1)
        plt.plot(dicct['x_sol_start'][0],dicct['x_sol_start'][1], "b*")
        # Plot the points where the bread leaves the table.
        plt.plot(dicct['table_traje'][0][index], dicct['table_traje'][1][index], "bo")
        # Plot the trajectroy up until this point
        plt.plot(dicct['table_traje'][0][:index+1], dicct['table_traje'][1][:index+1], label=r"$v$={}".format(v))
        # Draw the bread 
        draw_bread(dicct['table_traje'][:,0:index:100], np.pi/6.5 + dicct['y'][0:index:100,2], "g")

        x_p_free_fall = np.array([dicct['x_free_fall'],dicct['free_fall_traj']])

        spec_index = find_nearest(x_p_free_fall[1], table_hight)

        delta_T = T[-1]/N
        T_fall = np.linspace(0, spec_index * delta_T, spec_index)

        angle_free_fall = np.pi/6.5 + winkel(T_fall)

        draw_bread(x_p_free_fall[:,0:spec_index:100000], angle_free_fall[0:spec_index:100000], "r")


        # Plot the trajectroy after the bread has left the table
        plt.plot(dicct['x_free_fall'], dicct['free_fall_traj'])


        plt.ylim(-1.5, 0.1)
        plt.xlim(-0.1, 1.5)
        plt.plot(0.0, 0.0, "b*")
        plt.xlabel(r"$x$ [m]")
        plt.ylabel(r"$y$ [m]")
        plt.legend()


    return T_fall, winkel(T_fall)
N=1000000
ta, winekla = test(np.linspace(0,1,N))
tb, wineklb = test(np.linspace(0,10,N))
plt.show()

Note that the huge test function is just to play around with the different T.


EDIT: To illustrate the problem, two plots with different time samplings (N=1000,10000,...) and different time lengths (the first picture for T_end = 7s and the second one for T_end=10s). enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Do you get a discontinuity when you change the step? Does the result converge for vanishingly small steps? $\endgroup$ – Mathieu Bouville Apr 20 at 10:54
  • 3
    $\begingroup$ Give more information on the big picture: what you are trying to do, what results you get. Then move on to the code. Too many trees, too little forest. $\endgroup$ – Mathieu Bouville Apr 20 at 11:00
  • $\begingroup$ @MathieuBouville I edited in as much context as I could (a full mathematical description of the problem).. To answer your fist questions: No discontinutiy, but not sure what you meant with your second question... $\endgroup$ – Sito Apr 20 at 12:58
  • $\begingroup$ If you try steps of 1, 0.1, 0.01, 0.001, etc. do you eventually get a step-independent result? $\endgroup$ – Mathieu Bouville Apr 20 at 13:52
  • $\begingroup$ @MathieuBouville The results seem to converge to some value for decreasing step size, but changing the end time produces completely different resutls $\endgroup$ – Sito Apr 20 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.