双标量中遇到的溢出和无效值 - 非线性 PDE 求解

问题描述

我正在寻找一维非线性偏微分方程的有限差分解

u_t = u_xx + u(u_x)^2

代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import math

''' 
We explore three different numerical methods for solving the PDE,with solution u(x,t),u_t = u_xx + u(u_x)^2 

for (x,t) in (0,1) . (0,1/5)
u(x,0) = 40 * x^2 * (1 - x) / 3
u(0,t) = u(1,t) = 0

'''
M = 30
dx = 1 / M
r = 0.25
dt = r * dx**2
N = math.floor(0.2 / dt)

x = np.linspace(0,1,M + 1)
t = np.linspace(0,0.2,N + 1)

U = np.zeros((M + 1,N + 1))                      # Initial array for solution u(x,t)

U[:,0] = 40 * x**2 * (1 - x) / 3                 # Initial condition (: for the whole of that array)
U[0,:] = 0                                       # Boundary condition at x = 0
U[-1,:] = 0                                      # Boundary condition at x = 1 (-1 means end of the array)




'''
Explicit Scheme - Simple Forward Difference Scheme
'''

for q in range(0,N - 1):
    for p in range(0,M - 1):
        b = 1 / (1 - 2 * r)
        C = r * U[p,q] * (U[p + 1,q] - U[p,q])**2
        U[p,q + 1] = b * (U[p,q] + r * (U[p + 1,q + 1] + U[p - 1,q + 1]) - C)


T,X = np.meshgrid(t,x)

fig = plt.figure()
ax = fig.gca(projection='3d')

surf = ax.plot_surface(T,X,U)
#fig.colorbar(surf,shrink=0.5,aspect=5) # colour bar for reference

ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('u(x,t)')

plt.tight_layout()
plt.savefig('FDExplSol.png',bBox_inches='tight')
plt.show()

我使用的代码产生以下错误：

overflow encountered in double_scalars
  C = r * U[p,q])**2

invalid value encountered in double_scalars
  U[p,q + 1]) - C)

invalid value encountered in double_scalars
  C = r * U[p,q])**2

Z contains NaN values. This may result in rendering artifacts.
  surf = ax.plot_surface(T,U)

我已经查找了这些错误，并且我假设平方项生成的值对于 dtype 来说太小了。但是，当我尝试更改 dtype 以考虑更大范围的数字 (np.complex128) 时，我遇到了同样的错误。

结果图显然缺少大部分内容。所以，我的问题是，我该怎么办？

解决方法

离散化表达式不正确。

应该

for q in range(0,N - 1):
    for p in range(0,M - 1):
        U[p,q + 1] = r * (U[p + 1,q] - 2 * U[p,q] + U[p - 1,q]) + r * U[p,q] * (U[p + 1,q] - U[p,q])