A Simple VBA Tutorial: Euler solution of a mass spring damper system. This tutorial demonstrates the use on the Euler numerical method to solve the mass spring damper problem mx(t)”+cx(t)’+kx(t)=0.

1) The Theory :-

To solve using Euler

mx”(t)+cx'(t)+kx(t)=0.                (1)

it must be reduced to first order differential equations hence we define:-

X=x(t) and                            (2)

V=X'(t)                                 (3)

Then equation (1) above can be written as:-

mV'(t)+cV(t)+kX(t) = 0          (4)

hence:-

V'(t)=-(c/m)V(t)-(k/m)X(t)      (5)

X'(t)=V(t) from (3)                 (6)

We start at our initial condition which is with mass displaced by 1 and velocity 0 hence:-

X(t0) = 1

V(t0) = 0

To find the values at a small time increment (dt=0.01) after the mass has been released we calculate the derivatives V'(t) and X'(t) from equations 5 and 6.

with m=6,c=0.5,k=2

X'(t0) = 0

V'(t0) = -(0.5/6)*0-(2/6)*1 = -0.3333

The new values of X and V at time 0+dt are then estimated as

X(t0.01)    = X(t0)+dt*X'(t0)

= 1+0.01*0 = 1

V(t0.01)    = V(t0)+dt*V'(t0)

=0+0.01*-0.0333 = -0.003333

Then repeat this to find the derivatives at t0.01 and so on.

2) The code :-

Option Explicit

‘Global declarations

Dim M As Double

Dim C As Double

Dim K As Double

‘Euler Method for solving

‘mx”+cx’+kx=0

‘a simple mass spring damper problem

‘************************************************************

‘ Calculate the derivatives

‘************************************************************

Public Sub CalcDer(inX As Double, inV As Double, inDx As Double, inDv As Double)

inDx = inV

inDv = -(K / M) * inX – (C / M) * inV

End Sub

‘************************************************************

‘MAIN

‘************************************************************

Private Sub CommandButton1_Click()

Dim X As Double

Dim V As Double

Dim dt As Double

Dim T As Double

Dim i As Integer

Dim iROw As Integer

Dim dX As Double

Dim dV As Double

iROw = 1

Worksheets(1).Cells(iROw, 1) = “T”

Worksheets(1).Cells(iROw, 2) = “X”

Worksheets(1).Cells(iROw, 3) = “V”

iROw = 2

‘************************************************************

‘Problem setup

‘************************************************************

T = 0                                        ‘Start at time 0

dt = 0.01                                   ‘Time step

X = 1                                        ‘Initial displacement

V = 0                                        ‘Initial Veclocity

‘************************************************************

‘ Solution 1000 * dt

‘ This bit solves for the new X & V values at each new time

‘ step

‘************************************************************

For i = 0 To 10000

‘Report data to Excel

Worksheets(1).Cells(iROw, 1) = T

Worksheets(1).Cells(iROw, 2) = X

Worksheets(1).Cells(iROw, 3) = V

‘calculate the previous time step derivatives so we can estimate

‘the values at this time step

CalcDer X, V, dX, dV

‘calculate the new values at this time step

X = X + dt * dX

V = V + dt * dV

T = T + dt

iROw = iROw + 1

Next i

End Sub

The solution is very mach an approximation and for some values of m, c and k will be unstable, a more accurate solution is given by the Runge Kutta method.

The next tutorial 9 calculates the frequency response of a mass spring damper system.