El problema proviene del hecho de que los derivados no se calculan correctamente y no se consideran en el orden correcto.
Puede romper el algoritmo de la siguiente manera: –
theta (i + 1) = theta (i) + dtheta (i) * dt ………. (1)
d2theta (i) = [dtheta (i) -dtheta (i-1)] / dt …… .. (2)
Al organizar (1) y (2) puedes escribir
theta (i + 1) = theta (i) + [d2theta (i) * dt + dtheta (i-1)] * dt ……. (3)
donde d2theta (i) = f (theta (i), dtheta (i))
Entonces en tu código
Usted coloca la condición inicial como
theta1 = const1 theta2 = const2
dtheta1 = 0.d0 y dtheta2 = 0.d0 ya que están en reposo
luego estimar d2theta1 y d2theta2 utilizando los valores anteriores para obtener
El d2theta.
¡Luego pon todos estos valores en la ecuación (3) para obtener el theta (i + 1)!
Para el caso inicial, considere dtheta (i-1) = 0.d0 !!
En lo que respecta al método de Euler, debe ser muy pequeño
dt tener estabilidad numérica o de lo contrario la solución diverge.
Creo que debería ser al menos alrededor de ~ 0.1 * sqrt [(l1 + l2) / g] para numérico
estabilidad. Creo que alrededor de 0.001 debería ser lo suficientemente bueno.
Puede probar diferentes valores de dt en orden sistemático para tener una idea.
- Se suspende una bola de masa con una cuerda en el medio de un automóvil. ¿En qué dirección se mueve la pelota cuando el automóvil acelera a un ritmo constante?
- Si aceleras un objeto, ¿aumenta su masa? ¿Como y por qué?
- ¿Es posible hacer que el aire sea tan denso que podamos flotar sobre él?
- ¿Por qué los planetas no caen al sol, a pesar de girar a su alrededor? ¿Por qué giran en órbitas elípticas?
- Júpiter se compone de, yo diría, partículas densas, pero ¿qué tamaño y espacio tienen estas partículas individuales entre sí?