Python 中的欧拉方法

4

我正在尝试在 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是否合适......)

1
你的代码存在一些问题,但我想先看到完整的错误回溯信息,复制并粘贴在你的问题中,以及你如何调用“Euler”。你能否请提供这些信息来完善你的问题?谢谢。 - gboffi
2个回答

3
你尝试使用的公式不是欧拉方法,而是当n趋近于无穷大时e的精确值。参见维基百科
$n = \lim_{n\to\infty} (1 + \frac{1}{n})^n$

欧拉法 用于解决一阶微分方程。

以下是两个指南,展示了如何实现欧拉法来解决一个简单的测试函数: 初学者指南数值ODE指南

为了回答这篇文章的标题,而不是你正在问的问题,我使用欧拉法来解决通常的指数衰减:

$\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

注意:我不确定如何正确显示LaTeX。

1
你确定你不是在尝试实现牛顿迭代法吗?因为牛顿迭代法用于近似根。
如果你决定使用牛顿迭代法,这里提供了一个略微改变的代码版本,可以近似计算平方根 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]

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

此外,一个大问题是在Python中使用^而不是**来表示乘方运算,这是合法的,但是它们是完全不同的(按位)操作。


我肯定是指欧拉方法,但是... ** 绝对是一个问题。谢谢。 - newpythonuser

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接