python方程求根公式,用Python、Euler公式求解非线性方程组

正如DYZ所说，在第一次循环迭代中，您的计算是不正确的，因为j[-1]是{}的最后一个元素，您已经将其初始化为零。在

但是，您的代码浪费了大量的RAM。我假设您只需要包含T结果的数组，加上初始值，而不是每个步骤计算的结果。下面的代码通过一个双for循环来实现这一点。在这段代码中，我们并没有从Numpy中得到任何好处，所以我不必费心导入它。在

注意，Euler积分不是很精确，通常需要使用比更复杂的积分算法所要求的更小的步长。正如DYZ所提到的，在循环结束之前，使用当前步长，计算会发生偏差。在

下面是一个使用较小步长的代码修改版本。在T = 100

dt = 0.00001

steps = int(T / dt)

substeps = int(steps / T)

# Recalculate `dt` to compensate for possible truncation

# in the `steps` and `substeps` calculations

dt = 1.0 / substeps

print('steps, substeps, dt:', steps, substeps, dt)

a = 0.3025

C0 = 0.5

mu = 50

#dj = -mu*(J**3 - (C - C0)*J - F)

#dc = J + C*F + a*J**2

#df = J*F - C

# Initial condition

j = 0.1

c = -0.5

f = 0.1

jlst, clst, flst = [j], [c], [f]

for i in range(T):

for _ in range(substeps):

j1 = j + (-mu * (j**3 - (c - C0)*j - f))*dt

c1 = c + (j + c * f + (a * j)**2)*dt

f1 = f + (j * f - c)*dt

j, c, f = j1, c1, f1

jlst.append(j)

clst.append(c)

flst.append(f)

def round_seq(seq, places=6):

return [round(u, places) for u in seq]

print('j:', round_seq(jlst), end='\n\n')

print('c:', round_seq(clst), end='\n\n')

print('f:', round_seq(flst), end='\n\n')

输出

^{pr2}$

在我的旧2GHz机器上需要75秒。在

使用dt = 0.000005(在这台机器上大约需要2分钟)j，c，和{}的最终值分别是-0.524774、-0.529217、-0.674293，所以看起来我们开始收敛了。在

感谢LutzL指出dt可能需要调整，因为steps和{}计算中的取整。在