Python — реализация чехарды для нелинейной силы пружины, переполнение встречается в double_scalars ⇐ Python

[Anonymous]

Я пытаюсь справиться с тривиальным уравнением SHM, в котором нелинейный член альфа включен в силу пружины.
(https://i.sstatic.net/lrrDP.png)
На данный момент у меня есть следующий код для чехарды, которую я пытался решить, используя:
(https://i.sstatic.net/NNlrU.png)
(https://i.sstatic.net/85BKE.png)
Код:

Код: Выделить всё

import matplotlib.pyplot as plt
import numpy

k      = 1.0 # Spring constant
m      = 1.0 # Mass
cycles = 2.0 # No. of periods to integrate over
alpha  = 0.5 # Nonlinear Parameter
x0     = 1.0 # Initial displacement
v0     = 0.0 # Initial velocity

def leapfrog( steps ):
"""Solve the nonlinear oscillator equation for several
oscillation cycles"""
omega = (k/m)**0.5
delta = 2.0*cycles*numpy.pi/omega/steps
x     = numpy.empty( steps+1 )
v     = numpy.empty( steps+1 )
t     = numpy.empty( steps+1 )
t[0]  = 0.0
x[0]  = x0
v[0]  = v0 + 0.5*delta*x[0]*(omega**2*(1+alpha*x[0])*x[0])
for i in range(steps):
t[i+1] = t[i] + delta
x[i+1] = x[i] + delta*v[i]
v[i+1] = v[i] - delta*(omega**2*(1+alpha*x[i+1])*x[i+1])
return t, x, v

def manufactured_solution(t, i):
omega = numpy.sqrt(k/m)
M = alpha*omega**2*x0**2*numpy.cos(omega*t[i])**2
x_exact = x0*numpy.cos(omega*t[i]) + M/omega**2
return x_exact

def l2_error_norm( t , x ):
"""Calculate the L2 relative error norm."""
steps  = len( x ) - 1
omega  = (k/m)**0.5
l2_err = 0.0
l2     = 0.0
for i in range(steps):
x_exact = manufactured_solution(t, i)
l2_err += (x[i] - x_exact)**2.0
l2     += x[i]**2.0
return (l2_err/l2)**0.5

n = 14
steps = 8
delta    = numpy.empty( n )
l2_error = numpy.empty( n )

for i in range(0,n):
t, x, v     = leapfrog( steps )
delta[i]    = (k/m)**0.5*(t[1]-t[0])
l2_error[i] = l2_error_norm( t , x )
plt.plot( t , x )
# plt.show()
steps *= 2

# Switch to a new plotting window, and plot the L2 error norm,
# with guidelines for first, second, and third order accuracy.
plt.figure()
plt.loglog( delta , l2_error , 'o' )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**1.0 )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**2.0 )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**3.0 )
plt.show()

Если я попытаюсь решить эту проблему, возникнут следующие ошибки:

Код: Выделить всё

nhm.py:38: RuntimeWarning: overflow encountered in double_scalars
v[i+1] = v[i] - delta*(omega**2*(1+alpha*x[i+1])*x[i+1])
nhm.py:49: RuntimeWarning: overflow encountered in double_scalars
l2_err += (x[i] - x_exact)**2.0
nhm.py:50: RuntimeWarning: overflow encountered in double_scalars
l2     += x[i]**2.0
nhm.py:51: RuntimeWarning: invalid value encountered in
double_scalars return (l2_err/l2)**0.5

Любая полезная информация приветствуется.

Подробнее здесь: https://stackoverflow.com/questions/755 ... e-overflow