时间戳

首页 > 解决方案 > Numpy float128 没有给出正确答案

numpy - Numpy float128 没有给出正确答案

问题描述

我在 Python 中创建了一个微分方程求解器（Runge-Kutta 4 阶方法）。然后我决定通过将参数 mu 设置为 0 并查看它返回的数值解来检查它的结果。

问题是，我知道这个解决方案应该会产生稳定的振荡，但我得到了一个发散的解决方案。

代码如下所示。我尝试通过使用 numpy float128 数据类型来解决这个问题（浮点精度的舍入误差）。但是求解器一直给我错误的答案。

代码是：

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

def f(t,x,v):
  f = -k/m*x-mu/m*v
  return(f)

def g(t,x,v):
  g = v
  return(g)

def srunge4(t,x,v,dt):
  k1 = f(t,x,v)
  l1 = g(t,x,v)

  k2 = f(t+dt/2, x+k1*dt/2, v+l1*dt/2)
  l2 = g(t+dt/2, x+k1*dt/2, v+l1*dt/2)

  k3 = f(t+dt/2, x+k2*dt/2, v+l2*dt/2)
  l3 = g(t+dt/2, x+k2*dt/2, v+l2*dt/2)  

  k4 = f(t+dt/2, x+k3*dt, v+l3*dt)
  l4 = g(t+dt/2, x+k3*dt, v+l3*dt)

  v = v + dt/6*(k1+2*k2+2*k3+k4)
  x = x + dt/6*(l1+2*l2+2*l3+l4)
  t = t + dt
  return([t,x,v])

mu = np.float128(0.00); k = np.float128(0.1); m = np.float128(6)

x0 = np.float128(5); v0 = np.float128(-10)

t0 = np.float128(0); tf = np.float128(1000); dt = np.float128(0.05)

def sedol(t, x, v, tf, dt):
  sol = np.array([[t,x,v]], dtype='float128')
  while sol[-1][0]<=tf:
    t,x,v = srunge4(t,x,v,dt)
    sol=np.append(sol,np.float128([[t,x,v]]),axis=0)
  sol = pd.DataFrame(data=sol, columns=['t','x','v'])
  return(sol)

ft_runge = sedol(t0, x0, v0, tf, dt=0.1)

plt.close("all")
graf1 = plt.plot(ft_runge.iloc[:,0],ft_runge.iloc[:,1],'b')
plt.show()

我是否以错误的方式使用 numpy float128？

标签： numpyfloating-pointprecisionodesolver

解决方案

您正在混合和到和srunge4的关联。根据函数关联和最终求和，关联应该是和。在方法的各个阶段的更新中必须遵守这一点。klxv(v,f,k)(x,g,l)

在第 4 阶段，它应该t+dt在第一个参数中。然而，由于t没有在导数计算中使用，这个错误在这里没有任何后果。

此外，如果您在dt=0.1.

带有这些更正和一些进一步简化的代码是

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

mu = np.float128(0.00); k = np.float128(0.1); m = np.float128(6)
x0 = np.float128(5); v0 = np.float128(-10)
t0 = np.float128(0); tf = np.float128(1000); dt = np.float128(0.05)

def f(t,x,v): return -(k*x+mu*v)/m
def g(t,x,v): return v

def srunge4(t,x,v,dt): # should be skutta4, Wilhelm Kutta gave this method in 1901
  k1, l1 = (fg(t,x,v) for fg in (f,g)) 
  # here is the essential corrections, x->l, v->k
  k2, l2 = (fg(t+dt/2, x+l1*dt/2, v+k1*dt/2) for fg in (f,g)) 
  k3, l3 = (fg(t+dt/2, x+l2*dt/2, v+k2*dt/2) for fg in (f,g))
  k4, l4 = (fg(t+dt  , x+l3*dt  , v+k3*dt  ) for fg in (f,g))

  v = v + dt/6*(k1+2*k2+2*k3+k4)
  x = x + dt/6*(l1+2*l2+2*l3+l4)
  t = t + dt
  return([t,x,v])

def sedol(t, x, v, tf, dt):
  sol = [[t,x,v]]
  while t<=tf:
    t,x,v = srunge4(t,x,v,dt)
    sol.append([t,x,v])
  sol = pd.DataFrame(data=np.asarray(sol) , columns=['t','x','v'])
  return(sol)

ft_runge = sedol(t0, x0, v0, tf, dt=2*dt)

plt.close("all")
fig,ax = plt.subplots(1,3)
ft_runge.plot(x='t', y='x', ax=ax[0])
ft_runge.plot(x='t', y='v', ax=ax[1])
ft_runge.plot.scatter(x='x', y='v', s=1, ax=ax[2])
plt.show()

它产生预期的椭圆，而幅度没有视觉可识别的变化。

推荐阅读