Solución numérica aproximada de modelos matemáticos

El modelo anterior se resolverá mediante el método de Euler $$ \begin{equation*} \frac{d v}{d t} \approx \frac{\Delta v}{\Delta t} = \frac{v(t_{i+1}) - v(t_{i})}{t_{i+1} - t_{i}} \end{equation*} $$

reemplazando en la ecuación diferencial $$ \begin{equation*} \frac{v(t_{i+1}) - v(t_{i})}{t_{i+1} - t_{i}} = g - \frac{c v(t_{i})}{m} \end{equation*} $$

despejando \( v(t_{i+1}) \) $$ \begin{equation*} v(t_{i+1}) = v(t_{i}) + (g - \frac{c}{m} v(t_{i})) (t_{i+1} - t_{i}) \end{equation*} $$

m = 68.1
c = 12.5
g = 9.81
ti = 0
vi = 0

println("t", '\t', "v")
println(ti, '\t', vi)

for i = 1:20
    tf = ti + 1
    vf = vi + ((g - ((c / m) * vi)) * (tf  - ti))
    println(tf, '\t', vf)
    vi = vf
    ti = tf
end

    t	v
    0	0
    1	9.81
    2	17.819339207048458
    3	24.358535387839858
    4	29.697438583904496
    5	34.05636689082364
    6	37.615198225107115
    7	40.52079326455148
    8	42.89305588119034
    9	44.82988115997332
    10	46.41119519081522
    11	47.702253342280855
    12	48.75633312526895
    13	49.61693277187891
    14	50.31956625721685
    15	50.89322883849129
    16	51.361593589135325
    17	51.74398830478596
    18	52.05619309465638
    19	52.31109157214236
    20	52.51920251705015

Comparamos las soluciones:

using PyCall
using PyPlot

m = 68.1
c = 12.5
g = 9.81

ta = linspace(0,20,21)
va = zeros(21)

for i = 1:21
    va[i] = ((g * m) / c) * (1 - exp(-ta[i] * (c/m)))
end

tn = zeros(21)
vn = zeros(21)

for i = 1:20
    tn[i+1] = tn[i] + 1
    vn[i+1] = vn[i] + ((g - ((c / m) * vn[i])) * (tn[i+1]  - tn[i]))
end

plot(ta, va, color="blue", linewidth=2.0,label="Solución exacta")
plot(tn, vn, color="red", linewidth=2.0,label="Solución numérica")
xlabel(L"$t$")
ylabel(L"$v$")
legend(loc="upper left",fancybox="true")
title("Caída libre")
grid("on")