计算rk方法/无打印的函数

Question

import math 
import matplotlib
import numpy as np
from numpy import linspace

tmax=10.0
n=2000
G=4
D=-1
m=2
t=np.linspace (0,400,n+1)
phi=10
dphi=delta=phi_dot=np.linspace(0,400,n+1)

def f(delta_dot,t):
    return ((G)*(D*delta+m))

def iterate (func,phi,delta,tmax,n):
    dt=tmax/(n-1)
    t=0.0

    for i in range(n):
        phi,delta = func (phi,delta,t,dt)
        t += dt

    return phi

def rk_iter(phi,delta,t,dt):
    k1=f(t,phi)
    k2=f(t+dt*0.5,phi+k1*0.5*dt)
    k3=f(t+dt*0.5,phi*k2*0.5*dt)
    k4=f(t*dt,phi*k3*dt)

    delta +=dt*(k1+2*k2+2*k3+k4)/6
    k1=k2=k3=k4=delta=phi_dot
    phi += dt*(k1+2*k2+2*k3+k4)/6
    return phi,delta

runge_kutta = lambda delta, phi,tmax,n:iterate(rk_iter,delta,phi,tmax,n)

def plot_result (delta,phi,tmax,n):
 dt=tmax/(n-1)
 error_rk=[]
 r_rk=[]
 t=0.0

 phi=phi_rk=phi
 delta=delta_rk=delta

 for i in range(n):
     phi_rk,delta_rk=rk_iter(phi_rk,delta_rk,t,dt=tmax/(n-1))
     t+=dt
     _plot("error.png","Error","time t", "error e",error_rk)

def _plot(title,xlabel,ylabel,rk):
    import matplotlib.pyplot as plt

plt.title(title)
plt.ylabel(ylabel)
plt.xlabel(xlabel)

plt.plot(rk,"r--",label="Runge-Kutta")

plt.legend(loc=4)
plt.grid(True)
plt.plot(runge_kutta,t)


print "runge_kutta=", runge_kutta(phi,delta,tmax,n)
print "tmax=",t

我不知道如何让函数plt.show()工作。打开绘图窗口需要做什么？你知道吗

Answer 1

您尚未定义f；相反，f是由语句from matplotlib import *从matplotlib导入的：

In [10]: import matplotlib

In [11]: matplotlib.f
Out[11]: Forward: "a"

In [12]: matplotlib.f(1,1)
                                     -
TypeError                                 Traceback (most recent call last)
<ipython-input-12-5843007d7dbe> in <module>()
  > 1 matplotlib.f(1,1)

TypeError: __call__() takes at most 2 arguments (3 given)

如果在脚本中永远不要使用import语句的*形式，将来您将省去许多类似这样的麻烦。例如，不要使用from math import *，而是使用以下两种形式之一：

import math

# Now refer to math.sin, math.pi, etc.

或者，只显式导入要使用的名称：

from math import sin, pi

Answer 2

在任何时候都不能调用plot_result过程。即使您调用了它，也不会填充rk和error_rk列表。您只需使用runge_kutta调用的结果，。。。你知道吗

正如在另一篇重复文章中所评论的，您将系统方程定义为f(y,t)，但将其用作f(t,y)。你知道吗

delta的用法有些混乱，有时是积分变量，有时是rk4步骤更新。你知道吗

在rk4步骤中，有些乘法放错了位置，应该有加法。那条线呢

    k1=k2=k3=k4=delta=phi_dot

完全是无稽之谈，使先前的计算无效，并使rk4在下一步的更新毫无意义。你知道吗

删除math和linspace的导入，两者都没有在代码中使用。将plt的别名移到顶部，并将其与不必要的matplotlib导入合并。你知道吗