python - 尝试通过 Euler 和 Runge_Kutta 方法求解二阶 DE

Question

我正在尝试使用 Euler 和 Range-Kutta 方法解决 Spring 质量问题并比较绘图。我已经为 Euler 和 Runge-Kutta 编写了函数，但是在调用函数解决我的问题之后，我的绘图似乎没有显示任何数据。请帮我修复情节并检查我的代码是否有任何错误，谢谢

#function Euler
def euler ( y, t, dt, derivative):
    y_next = y + derivative(y, t) * dt
    return y_next

# function Runge-Kutta
# 2nd order Runge-Kutta method routine

def Runge_Kutta (y, time, dt, derivative):
    k0 = dt * derivative (y, time)
    k1 = dt * derivative (y + k0, time + dt)
    y_next = y + 0.5 * (k0 + k1)
    return y_next

这是我要解决的问题

[![""" A spring and mass  system. the coefficient of friction \mu is not negligible.generate a position vs. time plot for the motion of the mass, given an initial displacement x = 0.2m , spring constant k = 42 N/m , mass m =0.25 Kg, coefficient of friction \mu = 0.15 and initial velocity v = 0

F = -kx +/-mu mg """


from pylab import *
from Runge_Kutta_routine import Runge_Kutta
from eulerODE import euler

N = 500      #input ("How many number of steps to take?")
x0 = 0.2
v0 = 0.0
tau = 3.0     #input ("What is the total time of the simulation in seconds?")
dt = tau /float ( N-1)

k = 41.0      #input (" what is the spring constant?")
m = 0.25      #input ("what is the mass of the bob?")
gravity = 9.8
mu = 0.15     #input ("what is the coefficient of friction?")

""" we create a Nx2 array for storing the results of our calculations. Each 2- element row will be used for the state of the system at one instant, and each instant is separated by time dt. the first element in each row will denote position, the second would be velocity"""

y = zeros (\[N,2\])

y \[0,0\] = x0
y \[0,1\] = v0

def SpringMass (state, time):
    """ we break this second order DE into two first order DE introducing dx/ dt = v & dv/dt = kx/ m +/- mu g....

Note that the direction of the frictional force changes depending on the sign of the velocity, we handle this with an if statement."""


    g0 = state\[1\]
    if g0 > 0:
        g1 = -k/m * state \[0\] - gravity * mu
    else:
        g1 = -k/m * state \[0\] + gravity * mu

    return array (\[g0, g1\])

# Now we do the calculations
# loop only N-1 so that we don;t run into a problem addresssing y\[N+1\] on the last point

    for j in range (N-1):
        #y \[j+1\] = euler ( y\[j\] , 0, dt, SpringMass)
        y \[j+1\] = Runge_Kutta ( y\[j\], 0 , dt, SpringMass)

# Now we plot the result

time = linspace ( 0 , tau, N)
plot ( time, y\[:,0\], 'b-', label ='position')

xlabel('time')
ylabel('position')


show()][1]][1] 

Answer 1

看起来你的循环以 for j in range (N-1): 开始计算数组 y 是缩进的，所以 Python 认为这些行是函数 SpringMass。由于这些行位于 return 语句之后，因此它们永远不会被执行。

要更正此问题，请移动这些行，使 for 行没有缩进，而其他行只有四个空格的缩进。看看是否能解决您的问题。

请注意，您在此处编写的代码仍然无法运行。你在方括号前有多余的反斜杠，你在一个未命名的模块中编写你的 Euler 和 Runge_Kutta 函数，但主代码期望它们在两个不同的模块中，等等。你也有很多坏作风的例子。这些问题可能是您尚未得到任何(其他)答案的原因。帮自己一个忙，在此处发布和 clean up your style 之前清理代码.