使用 Euler 绘制解决方案 2nd ODE

Question

我已将运动方程（牛顿定律）用于简单的 spring 和质量场景，并将其纳入给定的第二个 ODE 方程 y" + (k/m)x = 0; y (0) = 3;y'(0) = 0.

使用欧拉法和精确解来解决问题，我已经能够运行并收到一些不错的结果。然而，当我执行结果图时，我得到了这条横跨我所追求的振荡结果的对角线。

Current plot output with diagonal line

任何人都可以帮助指出导致此问题的原因，我该如何解决？

我的代码：

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
from sympy.abc import x, i
import math


# Given is y" + (k/m)x = 0; y(0) = 3; y'(0) = 0

# Parameters
h = 0.01;  #Step Size
t = 50.0;  #Time(sec)
k = 1;     #Spring Stiffness
m = 1;     #Mass
x0 = 3;
v0 = 0;

# Exact Analytical Solution
x_exact = x0*cos(math.sqrt(k/m)*t);
v_exact = -x0*math.sqrt(k/m)*sin(math.sqrt(k/m)*t);

# Eulers Method
x = np.zeros( int( t/h ) );
v = np.zeros( int( t/h ) );
x[1] = x0;
v[1] = v0;
x_exact = np.zeros( int( t/h ) );
v_exact = np.zeros( int( t/h ) );
te      = np.zeros( int( t/h ) );
x_exact[1] = x0;
v_exact[1] = v0;


#print(len(x));

for i in range(1, int(t/h) - 1):    #MAIN LOOP
    x[i+1] = x[i] + h*v[i];
    v[i+1] = v[i] - h*k/m*x[i];
    te[i]  = i * h
    x_exact[i] = x0*cos(math.sqrt(k/m)* te[i]);
    v_exact[i] = -x0*math.sqrt(k/m)*sin(math.sqrt(k/m)* te[i]);
#    print(x_exact[i], '\t'*2, x[i]);

#plot
%config InlineBackend.figure_format = 'svg'
plt.plot(te, x_exact, te ,v_exact)
plt.title("DISPLACEMENT")
plt.xlabel("Time (s)")
plt.ylabel("Displacement (m)")
plt.grid(linewidth=0.3)

Answer 1

在一些细节上更直接的计算是

te = np.arange(0,t,h)
N = len(te)

w = (k/m)**0.5
x_exact = x0*np.cos(w*te);
v_exact = -x0*w*np.sin(w*te);

plt.plot(te, x_exact, te ,v_exact)

导致

请注意 python 中的数组从索引零开始，

x = np.empty(N)
v = np.empty(N)
x[0] = x0;
v[0] = v0;

for i in range(N - 1):    #MAIN LOOP
    x[i+1] = x[i] + h*v[i];
    v[i+1] = v[i] - h*k/m*x[i];

plt.plot(te, x, te ,v)

然后给出情节

具有预期的增加幅度。

使用 Euler 绘制解决方案 2nd ODE

Plotting solution 2nd ODE using Euler

ode

python-3.x

eulers-number

jupyter-notebook