2013-03-21 60 views
2

我是Python中的初级/中级。我已经编写了一个4th-order Runge-Kutta method (RK4)到Python。它基本上解决了摆锤,但这不是重点。我希望能够将函数f直接传递给RK4函数,即RK4(y_0,n,h)应该成为RK4(f,y_0,n) ,H)。这可以带来很大的好处,我可以将RK4用于描述其他系统的其他f函数,而不仅仅是这一个钟摆。如何将函数传递给Python中的函数?

我曾经玩过简单的功能到RK4,但我做错了什么。我如何在Python中做到这一点?

import numpy as np 

def RK4(y_0, n, h): 
    #4th order Runge-Kutta solver, takes as input 
    #initial value y_0, the number of steps n and stepsize h 
    #returns solution vector y and time vector t 
    #right now function f is defined below 

    t = np.linspace(0,n*h,n,endpoint = False) #create time vector t 
    y = np.zeros((n,len(y_0))) #create solution vector y 
    y[0] = y_0 #assign initial value to first position in y 
    for i in range(0,n-1): 
     #compute Runge-Kutta weights k_1 till k_4 
     k_1 = f(t[i],y[i]) 
     k_2 = f(t[i] + 0.5*h, y[i] + 0.5*h*k_1) 
     k_3 = f(t[i] + 0.5*h, y[i] + 0.5*h*k_2) 
     k_4 = f(t[i] + 0.5*h, y[i] + h*k_3) 
     #compute next y   
     y[i+1] = y[i] + h/6. * (k_1 + 2.*k_2 + 2.*k_3 + k_4) 
    return t,y 

def f(t,vec): 
    theta=vec[0] 
    omega = vec[1] 
    omegaDot = -np.sin(theta) - omega + np.cos(t) 
    result = np.array([omega,omegaDot])  
    return result 

test = np.array([0,0.5]) 
t,y = RK4(test,10,0.1) 

回答

4

Python函数也是对象。您可以通过他们周围像任何其他对象:

>>> def foo(): print 'Hello world!' 
... 
>>> foo 
<function foo at 0x10c4685f0> 
>>> foo() 
Hello world! 
>>> bar = foo 
>>> bar() 
Hello world! 

简单地传递一个函数作为一个额外的参数你RK4功能和使用,作为一个局部变量。

+0

感谢您的澄清。我来自MatLab(与一些C)。所以Python不停地惊讶于:-)。 – seb 2013-03-21 10:57:29

2

这很简单。更改RK4函数的定义,像这样:

def RK4(f, y_0, n, h): 

在这里,我已经添加了一个额外的参数,函数。

然后,当你调用RK4,传递函数:

t, y = RK4(f, test, 10, 0.1) 

而现在,当然,您也可以替换不同的功能,而无需重新编写集成代码。

Python中的函数只是另一种对象。你可以将它们传递给你,就像你做更多的平行对象一样。

+0

好吧,我把它放在里面,现在它可以工作。我是一个血腥的新人!无论如何,将它们全部输入并张贴在这里是一个很好的练习。 – seb 2013-03-21 10:56:01

+0

是的,很多时候写一个简洁的例子有助于你的理解。顺便说一下,这是一个写得很好的问题,做得很好! – 2013-03-21 10:57:33

3

您可以传递一个函数的功能在Python,正如你所期望的:

def call_function(f): 
    f() 

def my_function(): 
    print "OK" 

call_function(my_function) # Prints OK 

也许你应该张贴发生故障的代码?