老师讲交通流理论留的作业，说要matlab做，我试了一下python，发现matlab在三维方面有很多函数还是比python好用。

![Uploading 01EC438C17A2DDDF44AB9CED0859B817_952446.jpg . . .]
![Uploading 52311C47E13ECBB2A6FB57393CABA79F_956948.jpg . . .]
![Uploading 9C1C388EAFB639D6A1E949179612704C_959808.jpg . . .]
![Uploading 0A0F932DF35A40D79BFADEADDC6A698C_961345.jpg . . .]

python代码

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

Uf=27.8#速度
Kjam=0.035#Kj阻塞密度
T=120#观测周期
t=0.3#时间间隔
L=2000#路段长度
x=10

K=np.zeros(shape=(200,400))
Z=[]
def k0j():
    for j in range(0,200):
        Z.append(j*10*(2000-j*10)/40000000)
    return Z
k0j()
K[:,0]=Z

for n in range(0,399):
    for j in range(0,200):
        if j==0:
            K[j,n+1]=1/2*(K[j+1,n]+0)-0.3/(2*10)*Uf*(K[j+1,n]*(1-K[j+1,n]/Kjam)-0)
        if j==199:
            K[j, n + 1] = 1 / 2 * (0 + K[j - 1, n]) - 0.3 / (2 * 10) * Uf * (0 - K[j - 1, n] * (1 - K[j - 1, n] / Kjam))
        else:
            K[j, n + 1] = 1 / 2 * (K[j + 1, n] + K[j - 1, n]) - 0.3 / (2 * 10) * Uf * (K[j + 1, n] * (1 - K[j + 1, n] / Kjam) - K[j - 1, n] * (1 - K[j - 1, n] / Kjam))
print(K)


X = [x for x in range(0,200) for y in range(0,400)]
Y = [x for x in range(0,400)]*200

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
Z=K.reshape(1,80000)


ax.scatter(X, Y, Z)

plt.xlabel('t/s')
plt.ylabel('x/m')

plt.show()

matlab代码

Uf=27.8;%速度;
Kjam=0.035;%Kj阻塞密度;
T=120;%观测周期
t=0.3;%时间间隔;
L=2000;%路段长度;
x=10;

K=zeros(200,400);
for j=1:200
    K(j,1)=j*10*(2000-j*10)/40000000;
end;

for n =1:399
    for j =1:200
        if j==1
            K(j,n+1)=1/2*(K(j+1,n)+0)-0.3/(2*10)*Uf*(K(j+1,n)*(1-K(j+1,n)/Kjam)-0);
        elseif j==200
            K(j, n + 1) = 1 / 2 * (0 + K(j - 1, n)) - 0.3 / (2 * 10) * Uf * (0 - K(j - 1, n) * (1 - K(j - 1, n) / Kjam));
        else
            K(j, n + 1) = 1 / 2 * (K(j + 1, n) + K(j - 1, n)) - 0.3 / (2 * 10) * Uf * (K(j + 1, n) * (1 - K(j + 1, n) / Kjam) - K(j - 1, n) * (1 - K(j - 1, n) / Kjam));
        end;
    end;
end;
surf(K)