使用Euler绘制解决方案第二ODE

问题描述

我已经将运动方程（牛顿定律）用于简单的弹簧和质量场景，并将其结合到给定的第二ODE方程y“ +（k / m）x = 0; y（0）= 3; y' （0）= 0。

使用Euler方法和精确的解决方案来解决问题，我已经能够运行并收到一些不错的结果。但是，当我执行结果绘图时，我得到的对角线跨越了我所跟踪的振荡结果。

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)

解决方法

更详细地讲，直接计算是

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,v)

然后给出情节

具有预期的幅度增加。