python中的Euler方法

Question

我试图实现euler's method来近似python中e的值。这就是我目前所拥有的：

def Euler(f, t0, y0, h, N):
    t = t0 + arange(N+1)*h
    y = zeros(N+1)
    y[0] = y0
    for n in range(N):
        y[n+1] = y[n] + h*f(t[n], y[n])
        f = (1+(1/N))^N
    return y

但是，当我尝试调用该函数时，得到错误“ValueError:shape<；=0”。我怀疑这和我怎么定义f有关？我试着在调用euler时直接输入f，但给出了与未定义变量相关的错误。我还尝试将f定义为它自己的函数，这给了我一个0的除法错误。

def f(N):
    for n in range(N): 
        return (1+(1/n))^n

（不确定N是否是此处使用的适当变量…）

Answer 1

你试图使用的公式不是欧拉方法，而是当n接近无穷大时e的精确值wiki

$n = \lim_{n\to\infty} (1 + \frac{1}{n})^n$

---

Euler's method用于求解一阶微分方程。

这里有两个指南，说明如何实现Euler方法来解决一个简单的测试函数：beginner's guide和numerical ODE guide。

为了回答这篇文章的标题，而不是你所问的问题，我使用了欧拉方法来解决通常的指数衰减：

$\frac{dN}{dt} = -\lambda N$

有解决办法

$N(t) = N_0 e^{-\lambda t}$

代码：

import numpy as np
import matplotlib.pyplot as plt
from __future__ import division

# Concentration over time
N = lambda t: N0 * np.exp(-k * t)
# dN/dt
def dx_dt(x):
    return -k * x

k = .5
h = 0.001
N0 = 100.

t = np.arange(0, 10, h)
y = np.zeros(len(t))

y[0] = N0
for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + dx_dt(y[i-1]) * h

max_error = abs(y-N(t)).max()
print 'Max difference between the exact solution and Euler's approximation with step size h=0.001:'

print '{0:.15}'.format(max_error)

输出：

Max difference between the exact solution and Euler's approximation with step size h=0.001:
0.00919890254720457

---

注：我不知道如何使乳胶显示正确。

Answer 2

你确定你没有试图实现牛顿法吗？因为牛顿法是用来逼近根的。

如果你决定采用牛顿方法，这里有一个稍微改变的代码版本，它近似于2的平方根。你可以用函数及其导数来改变f(x)和fp(x)，你可以用它来近似你想要的东西。

import numpy as np

def f(x):
    return x**2 - 2


def fp(x):
    return 2*x

def Newton(f, y0, N):
    y = np.zeros(N+1)
    y[0] = y0
    for n in range(N):
        y[n+1] = y[n] - f(y[n])/fp(y[n])
    return y

print Newton(f, 1, 10)

给予

[ 1. 1.5 1.41666667 1.41421569 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356 1.41421356]

这是初始值和前十次迭代的平方根。

除此之外，一个大问题是使用^而不是**作为powers，这在python中是合法的但完全不同（按位）的操作。