Math Forum / Mathematics / Numerical Analysis / October 2008

Hi all,
It is a numeric calculus problem , not a tool problem.
You can solve your differenzial equation with the tool that you like
(Excel, Matlab, Scilab, Octave, ecc)

With Xnumbers you can solve easily your ODE with many algorithms:
(Runge-Kutta 4, Implicit A-stable algorithm, Runge-Kutte-Feldberg 4.5,
Exponential, Finite differences, etc.).
For using these feauture, please, consult the ODE_Tutoria pdf
at pag 38 it is showing how to solve just your 2nd order equation

http://digilander.libero.it/foxes/diffequ/ODE_tutorial.pdf

The problem is the integration step that depends on the parameter m,
k, b, g (that you have not mentioned)
I have try to solve the ODE using the function ODE_RK4 (Runge-Kutta 4
order) with this parameters: m = 2 kg (mass), k= 10 (elastic
constant), b = 3 (damping factor), g = 9,8 (acc. gravity) obtaining a
good result, with a step dt = 0.05 and 200 steps.
Of course I have also implentented you code in VBA Excel
I have only correct the assignement operator because you use two "="
for each line
and this is accepted in VB but has anothe meaning: it return TRUE or
FALSE, not the correct numerica value.

the correrect implementation is the following
(I don't know if it is allowed to attach file. If you want I can send
you the xls file with the macro at your private e-mail)
Hope this help you.
Good calculation.

Leonardo Volpi
Foxes Team, Italy
----------------------------
Code implementing Runge Kutta 4
<<<<...

Sub rk4_test()

'declaration
Dim imax, i, dt, t_i, x_i, v_i, v_i_1, x_i_1
Dim vk0, vk1, vk2, vk3, xk0, xk1, xk2, xk3

'parameter setting
k = 10
m = 1
b = 3
g = 9.8

dt = 0.05 'integration step
imax = 200 'max step
i = 0

'initial conditions
t_i = 0
x_i = 0
v_i = 0

'RK4 starts
Do
vk0 = f_v(x_i, v_i, t_i)
vk1 = f_v(x_i + xk0 * (dt / 2), v_i + vk0 * (dt / 2), t_i + dt / 2)
vk2 = f_v(x_i + xk1 * (dt / 2), v_i + vk1 * (dt / 2), t_i + dt / 2)
vk3 = f_v(x_i + xk2 * dt, v_i + vk2 * dt, t_i + dt)

xk0 = f_x(x_i, v_i, t_i)
xk1 = f_x(x_i + xk0 * (dt / 2), v_i + vk0 * (dt / 2), t_i + dt / 2)
xk2 = f_x(x_i + xk1 * (dt / 2), v_i + vk1 * (dt / 2), t_i + dt / 2)
xk3 = f_x(x_i + xk2 * dt, v_i + vk2 * dt, t_i + dt)

v_i_1 = v_i + (dt / 6) * (vk0 + 2 * vk1 + 2 * vk2 + vk3)
x_i_1 = x_i + (dt / 6) * (xk0 + 2 * xk1 + 2 * xk2 + xk3)

debug.print i, t_i, x_i, v_i

i = i + 1
v_i = v_i_1
x_i = x_i_1
t_i = t_i + dt

Loop Until i > imax

exit sub

Private Function f_v(x_i, v_i, t_i)
f_v = -(k / m) * x_i - (b / m) * v_i + g
End Function

Private Function f_x(x_i, v_i, t_i)
f_x = v_i
End Function

<<<<...