0

I am working on using the forward difference scheme for numerically solving the diffusion function in one dimension. My final plot of the solution should be a surface where the solution u(x,t) is plotted over a grid of x and t values. I have the problem solved, but I can't get the data to be plotted with the grid representation.

I can think of 2 ways to fix this:

1.) My x and t arrays should be one dimensional, but my u array should be a 2D array. Ultimately, I want a square matrix for u, but I am having a hard time coding that. Currently I have a 1D array for u. Here is the code where u is populated.

u   = zeros(Nx+1)           # unknown u at new time level
u_1 = zeros(Nx+1)           # u at the previous time level
# Set initial condition u(x,0) = I(x)
for i in range(0, Nx+1):
#set initial u's to I(xi)
    u_1[i] = 25-x[i]**2
for n in range(0, Nt):
# Compute u at inner mesh points
    for i in range(1, Nx):
        u[i] = u_1[i] + F*(u_1[i-1] - 2*u_1[i] + u_1[i+1])

2.) The above code returns a 1D array for u, is there a way to plot a 3D surface with 3 1D arrays for x,y,z?

4

1 回答 1

0

嗯...,有很多你没有提供的信息。例如,您说您想要 ax,y,z 绘图,但没有说明 x、y 和 z 在绘图上下文中应该是什么。z 通常也是 z(x,y)。

以下配方假设 atx, 和u(t,x)作为要放入表面的变量。我想这不完全是您的想法,但它应该适合您的锻炼:

编辑:你的代码(在这个配方的函数computeU中)也有一个循环Nt,似乎没有做任何事情。出于本示例的目的,我已将其删除。

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

def computeU(Nx,x,F,Nt):
    u   = np.zeros(Nx+1)           # unknown u at new time level
    u_1 = np.zeros(Nx+1)           # u at the previous time level
    # Set initial condition u(x,0) = I(x)
    for i in range(0, Nx+1):
    #set initial u's to I(xi)
        u_1[i] = 25-x[i]**2
    #for n in range(0, Nt): # I'm not sure what this is doing. It has no effect.
    # Compute u at inner mesh points
    for i in range(1, Nx):
        u[i] = u_1[i] + F*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
    return np.hstack((u[:,np.newaxis],u_1[:,np.newaxis]))

Nx = 10
F  = 3
Nt = 5
x  = np.arange(11)
t  = np.arange(2)

X,Y = np.meshgrid(t,x)
Z = computeU(Nx,x,F,Nt)
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)
plt.show()

注意我是如何meshgrid构建新t的 , x(从一维数组)映射到你的U数组堆栈(将具有相同的形状X, Y- new t, x)。结果是这样的:

由多个 1D 阵列构建的表面

于 2016-05-05T20:10:57.410 回答