ode - time and frequency domain ODE solvers =========================================== This and other notebooks are available here: https://github.com/twmacro/pyyeti/tree/master/docs/tutorials. Primarily, the `pyyeti.ode <../modules/ode.html>`__ module provides tools for solving 2nd order matrix equations of motion in both the time and frequency domains. This notebook demonstrates solving time-domain equations of motion. Notes: - Some features depend on the equations being in modal space (particularly important where there are distinctions between the rigid-body modes and the elastic modes). - The time-domain solvers all depend on constant time step. First, do some imports: .. code:: ipython3 import numpy as np import matplotlib.pyplot as plt from pyyeti import ode Some settings specifically for the jupyter notebook. .. code:: ipython3 %matplotlib inline plt.rcParams['figure.figsize'] = [6.4, 4.8] plt.rcParams['figure.dpi'] = 150. Setup a simple system ^^^^^^^^^^^^^^^^^^^^^ This system is not fixed to ground so it does have a rigid-body mode. This system is from the `frclim.ntfl <../modules/generated/pyyeti.frclim.ntfl.html#pyyeti.frclim.ntfl>`__ example. :: |--> x1 |--> x2 |--> x3 |--> x4 | | | | |----| k1 |----| k2 |----| k3 |----| Fe | |--\/\/\--| |--\/\/\--| |--\/\/\--| | ====>| 10 | | 30 | | 3 | | 2 | | |---| |---| |---| |---| |---| |---| | |----| c1 |----| c2 |----| c3 |----| Define the mass, damping and stiffness matrices: .. code:: ipython3 M1 = 10. M2 = 30. M3 = 3. M4 = 2. c1 = 15. c2 = 15. c3 = 15. k1 = 45000. k2 = 25000. k3 = 10000. MASS = np.array([[M1, 0, 0, 0], [0, M2, 0, 0], [0, 0, M3, 0], [0, 0, 0, M4]]) DAMP = np.array([[c1, -c1, 0, 0], [-c1, c1+c2, -c2, 0], [0, -c2, c2+c3, -c3], [0, 0, -c3, c3]]) STIF = np.array([[k1, -k1, 0, 0], [-k1, k1+k2, -k2, 0], [0, -k2, k2+k3, -k3], [0, 0, -k3, k3]]) Define a time vector and a forcing function: .. code:: ipython3 h = .001 t = np.arange(0, 2, h) F = np.zeros((4, len(t))) F[0, :200] = np.arange(200) F[0, 200:400] = np.arange(200)[::-1] F[0, 400:600] = np.arange(0, -200, -1) F[0, 600:800] = np.arange(0, -200, -1)[::-1] plt.plot(t, F[0]) plt.xlabel('Time (s)') plt.ylabel('Force') plt.grid(True) .. image:: temp_ode_files/temp_ode_11_0.png The `pyyeti.ode <../modules/ode.html>`__ module was originally designed only to solve systems that were in modal space via the `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ solver. Such equations are typically uncoupled, unless the damping is coupled. This is an exact solver assuming piece-wise linear forces. For the current problem, which is not in modal space, the solver now accepts the ``pre_eig=True`` option which will internally transform the problem to modal space. You can also choose the `pyyeti.ode.SolveExp2 <../modules/generated/pyyeti.ode.SolveExp2.html#pyyeti.ode.SolveExp2>`__ solver for the current problem. This solver is based on the matrix exponential and is more general than `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__. For uncoupled equations however `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ is likely significantly faster. Note, even for the `pyyeti.ode.SolveExp2 <../modules/generated/pyyeti.ode.SolveExp2.html#pyyeti.ode.SolveExp2>`__ solver, the equations will need to be put in modal space (via the ``pre_eig=True``) if you wish to use static initial conditions. The `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ solver ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Since this system is not in modal space, we’ll use the ``pre_eig=True``. .. code:: ipython3 ts = ode.SolveUnc(MASS, DAMP, STIF, h, pre_eig=True) At this point, many pre-calculations are done and the ``ts`` object is ready to be called (repeatedly if necessary) to solve the equations of motion with arbitrary forces and initial conditions. The time domain solver is called via the method `tsolve <../modules/generated/pyyeti.ode.SolveUnc.tsolve.html#pyyeti.ode.SolveUnc.tsolve>`__: .. code:: ipython3 sol = ts.tsolve(F) ``sol`` is a SimpleNamespace with the members: - .a - acceleration - .v - velocity - .d - displacement - .t - time vector - .h - time step Plot the acceleration responses (assume metric units): .. code:: ipython3 for j in range(sol.a.shape[0]): plt.plot(sol.t, sol.a[j], label='Acce {:d}'.format(j+1)) plt.legend(loc='best') plt.xlabel('Time (s)') plt.ylabel(r'Acceleration (m/$sec^2$)') plt.grid(True) .. image:: temp_ode_files/temp_ode_18_0.png Plot the velocities: .. code:: ipython3 for j in range(sol.v.shape[0]): plt.plot(sol.t, sol.v[j], label='Velocity {:d}'.format(j+1)) plt.legend(loc='best') plt.xlabel('Time (s)') plt.ylabel('Velocity (m/sec)') plt.grid(True) .. image:: temp_ode_files/temp_ode_20_0.png Plot the displacements: .. code:: ipython3 for j in range(sol.d.shape[0]): plt.plot(sol.t, sol.d[j], label='Displacement {:d}'.format(j+1)) plt.legend(loc='best') plt.xlabel('Time (s)') plt.ylabel('Displacement (m)') plt.grid(True) .. image:: temp_ode_files/temp_ode_22_0.png -------------- The `pyyeti.ode.SolveExp2 <../modules/generated/pyyeti.ode.SolveExp2.html#pyyeti.ode.SolveExp2>`__ solver ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For demonstration, we’ll solve the same system using the `pyyeti.ode.SolveExp2 <../modules/generated/pyyeti.ode.SolveExp2.html#pyyeti.ode.SolveExp2>`__ solver. Since we’re not using static initial conditions, we don’t need to use the ``pre_eig`` option. In this case, not using it would be more efficient since it would do less work. For demonstration, we’ll use both: .. code:: ipython3 ts1 = ode.SolveExp2(MASS, DAMP, STIF, h) ts2 = ode.SolveExp2(MASS, DAMP, STIF, h, pre_eig=True) Solve with both solvers and demonstrate they give the same results: .. code:: ipython3 sol1 = ts1.tsolve(F) sol2 = ts2.tsolve(F) assert np.allclose(sol1.a, sol2.a) assert np.allclose(sol1.v, sol2.v) assert np.allclose(sol1.d, sol2.d) Show the results are also the same as the `pyyeti.ode.SolveUnc <../modules/generated/pyyeti.ode.SolveUnc.html#pyyeti.ode.SolveUnc>`__ solver: .. code:: ipython3 assert np.allclose(sol1.a, sol.a) assert np.allclose(sol1.v, sol.v) assert np.allclose(sol1.d, sol.d)