352 - Topic 10 - Automobile Suspension Model -- Sage

# Solving a second-order linear DE in general using Sage var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

Traceback (click to the left of this block for traceback)
...
NotImplementedError: Maxima was unable to solve this ODE. Consider to
set option contrib_ode to True.

Traceback (most recent call last):    soln=desolve(DE,y,ivar=t)
  File "", line 1, in <module>
    
  File "/tmp/tmp3vEZ7E/___code___.py", line 8, in <module>
    soln=desolve(DE,y,ivar=t)
  File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/calculus/desolvers.py", line 597, in desolve
    raise NotImplementedError("Maxima was unable to solve this ODE. Consider to set option contrib_ode to True.")
NotImplementedError: Maxima was unable to solve this ODE. Consider to set option contrib_ode to True.

# Solving a second-order linear DE in general using Sage # Let's try this with specific values for a1 and a0. var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) a1 = t a0 = 1 DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} i \, \sqrt{2} \sqrt{\pi} K_{1} \operatorname{erf}\left(\frac{1}{2} i \, \sqrt{2} t\right) e^{\left(-\frac{1}{2} \, t^{2}\right)} + K_{2} e^{\left(-\frac{1}{2} \, t^{2}\right)}$

# Solving a second-order linear DE in general using Sage # Let's try this with specific constant values for a1 and a0. var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) a1 = 2 a0 = 2 DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln.expand())

$\newcommand{\Bold}[1]{\mathbf{#1}}K_{2} \cos\left(t\right) e^{\left(-t\right)} + K_{1} e^{\left(-t\right)} \sin\left(t\right)$

# Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) a1 = 1 a0 = 1 DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(K_{2} \cos\left(\frac{1}{2} \, \sqrt{3} t\right) + K_{1} \sin\left(\frac{1}{2} \, \sqrt{3} t\right)\right)} e^{\left(-\frac{1}{2} \, t\right)}$

# Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) a1 = 3 a0 = 1 DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

$\newcommand{\Bold}[1]{\mathbf{#1}}K_{2} e^{\left(-\frac{1}{2} \, t {\left(\sqrt{5} + 3\right)}\right)} + K_{1} e^{\left(\frac{1}{2} \, t {\left(\sqrt{5} - 3\right)}\right)}$

Automobile Suspension (Harmonic Oscillator)

Assume a wheel with mass m on a stationary automobile is suspended and the equilibrium position presume the position of the wheel is at y = 0.

At the equilibrium, the spring is not moving and so the force holding the wheel back is equal to the force of gravity pulling the wheel down. If $\lambda$ is the distance that gravity stretches the spring (and moves the wheel down) from it's natural position, then Hooke's Law suggests that the force exerted by the spring is proportional to this length or

$mg = k\lambda$

where k is the Hooke's Law proportionality constant.

Suppose now that the wheel is moved from this equilibrium position and that y is the distance the wheel is from equilibrium where positive y-values are oriented downward (so that y>0 means the spring is stretched and y<0 means the spring is compressed). Hooke's law again says that an additional force will be exerted proportional to this distance y.So, the net resulting force on the spring when the wheel is at location y is

gravity + wheel weight on spring + spring's pulling/pushing = $mg - k(y+\lambda)$

where we note that y>0 means the spring is pulling back up (in the negative direction) and y<0 means the spring is pushing down (in the positive direction).

Therefore the net force on the wheel by Newton's Second Law is again $F = ma = m y''$ and so we get the differential equation

$m y'' = mg - k(y+\lambda) = k\lambda - ky - k\lambda $

$my'' + ky = 0$

which is a second order constant coefficient homogeneous DE.

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t,c,k') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) assume(k>0) assume(c>0) a1 = 0 a0 = k DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

$\newcommand{\Bold}[1]{\mathbf{#1}}K_{2} \cos\left(\sqrt{k} t\right) + K_{1} \sin\left(\sqrt{k} t\right)$

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) assume(k>0) @interact def _(a0 = slider(0.1,2,0.1), y0 = input_box(default=1, width = 10), # initial starting position yp0 = input_box(default=0, width = 10)): # initial starting velocity DE = diff(y,t,2) + a0*y==0 t0 = 0 # startign time soln=desolve(DE,y,ivar=t, ics = [t0,y0,yp0]) show(soln) plot(soln,t,0,10*pi).show()

Click to the left again to hide and once more to show the dynamic interactive window

Notice, the solution to the DE above is unsatisfactory for safe travel since once a car hits a bump that it continues to bounce forever. Not a safe ride. So, let's add a shock absorber that will affect the motion of the wheel only while the wheel is moving. Since the first derivative models velocity (how fast a thing moves) then perhaps a good model would be for a shock absorber to exert a force that is negatively proportional to the speed of any movement. That is a force that looks like $ -c y' $. Adding this into our net force equation from before gives

$$ m y'' = mg - k(y+\lambda)- cy' = k\lambda - ky - k\lambda - cy' $$

$ m y'' + cy' + ky = 0 $

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t,c,k') y = function('y')(t) a1 = function('P')(t) a0 = function('Q')(t) assume(k>0) assume(c>0) a1 = c a0 = k DE = diff(y,t,2) + a1*diff(y,t) + a0*y==0 soln=desolve(DE,y,ivar=t) show(soln)

Traceback (click to the left of this block for traceback)
...
Is 4*k-c^2 positive, negative or zero?

Traceback (most recent call last):    a0 = function('Q')(t)
  File "", line 1, in <module>
    
  File "/tmp/tmp812o80/___code___.py", line 14, in <module>
    soln=desolve(DE,y,ivar=t)
  File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/calculus/desolvers.py", line 584, in desolve
    soln = P(cmd)
  File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/interfaces/interface.py", line 288, in __call__
    return cls(self, x, name=name)
  File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/interfaces/interface.py", line 703, in __init__
    raise TypeError(x)
TypeError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(4*k-c^2>0)', see `assume?` for more details)
Is 4*k-c^2 positive, negative or zero?

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t,c,k') y = function('y')(t) assume(k>0) assume(c>0) assume(4*k-c^2>0) DE = diff(y,t,2) + c*diff(y,t) + k*y==0 soln=desolve(DE,y,ivar=t) show(soln)

$\newcommand{\Bold}[1]{\mathbf{#1}}{\left(K_{2} \cos\left(\frac{1}{2} \, \sqrt{-c^{2} + 4 \, k} t\right) + K_{1} \sin\left(\frac{1}{2} \, \sqrt{-c^{2} + 4 \, k} t\right)\right)} e^{\left(-\frac{1}{2} \, c t\right)}$

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t,c,k') y = function('y')(t) assume(k>0) assume(c>0) assume(4*k-c^2<0) DE = diff(y,t,2) + c*diff(y,t) + k*y==0 soln=desolve(DE,y,ivar=t) show(soln)

Traceback (click to the left of this block for traceback)
...
ValueError: Assumption is inconsistent

Traceback (most recent call last):    assume(k>0)
  File "", line 1, in <module>
    
  File "/tmp/tmpmyN0Om/___code___.py", line 10, in <module>
    assume(_sage_const_4 *k-c**_sage_const_2 <_sage_const_0 )
  File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/symbolic/assumptions.py", line 550, in assume
    x.assume()
  File "sage/symbolic/expression.pyx", line 1837, in sage.symbolic.expression.Expression.assume (build/cythonized/sage/symbolic/expression.cpp:13610)
ValueError: Assumption is inconsistent

Notice the bifuraction when you pass across the point where $c^2-4k=0$.

# Solving basic Harmonic Oscillator # Solving a second-order linear DE in general using Sage # Let's try this with specific co values for a1 and a0. var('t') y = function('y')(t) assume(k>0) @interact def _(k = slider(0.1,2,0.1), # spring stiffness c = slider(0.1,2,0.1), # shock stiffness y0 = input_box(default=1, width = 10), # initial starting position yp0 = input_box(default=0, width = 10)): # initial starting velocity pretty_print(c^2-4*k) DE = diff(y,t,2) + c*diff(y,t) + k*y==0 t0 = 0 # startign time soln=desolve(DE,y,ivar=t, ics = [t0,y0,yp0]) show(soln) plot(soln,t,0,20*pi).show()

Click to the left again to hide and once more to show the dynamic interactive window

solve(x^2+c*x+k==0,x)

[x == -1/2*c - 1/2*sqrt(c^2 - 4*k), x == -1/2*c + 1/2*sqrt(c^2 - 4*k)]

[x == -1/2*c - 1/2*sqrt(c^2 - 4*k), x == -1/2*c + 1/2*sqrt(c^2 - 4*k)]

352 - Topic 10 - Automobile Suspension Model

10 hours ago by Professor352

Click to the left again to hide and once more to show the dynamic interactive window

Click to the left again to hide and once more to show the dynamic interactive window