我试图用欧拉和Range-Kutta方法来解决弹簧质量问题，并比较这些图。我已经为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]
最佳答案
看起来以for j in range (N-1):开始的计算数组y的循环是缩进的，因此Python将这些行视为函数SpringMass的一部分。因为这些行在return语句之后，所以它们永远不会被执行。
要更正此问题，请移动这些行，使for行没有缩进，而其他行只有四个缩进空间。看看能不能解决你的问题。
请注意，此处编写的代码仍然无法工作。在方括号前有多余的反斜杠，在一个未命名的模块中编写Euler和Runge_Kutta函数，但主代码希望它们位于两个不同的模块中，等等，还有许多样式不好的示例。这些问题可能是你还没有得到任何（其他）答案的原因。帮你自己一个忙，在这里和clean up your style发布之前清理你的代码。
关于python - 尝试通过Euler和Runge_Kutta方法求解二阶DE，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/37408700/