本人编程小白。最近在把之前写的解偏微分的代码重写成taichi，发现了两个问题：
第一，同样的公式数据计算结果与之前写的python代码的结果，和解析解相差很大。我演算了初始5次iteration中的微分du/dt结果正常，之后的微分值就和解析解差很多了。我不太清楚这是因为我的公式的问题还是其他原因。
第二，我发现同一个解微分方程的程序，每次运算出的结果都有不同，偶尔会出现完全错误的结果。这种现象正常吗？
最后附上代码，请带哥们不吝赐教！
import numpy as np
import taichi as ti
import time
ti.init(arch=ti.cpu, default_fp=ti.f32, dynamic_index=True)
@ ti.kernel
def inti_main():
t = 0. # end time
tmax = 0.01 # start time
dt = ti.cast(0.000001,dtype=ti.f32) # delta time
dx = ti.cast(2/32, dtype=ti.f32) # dx is prepared but unused in this case
leng = ti.cast(tmax/dt, dtype=ti.i32) # length of iteration---tmax/dt, better define manually
# 'leng' should be change along with the shape of field
u = init_con
for c in range(leng):
u_new, t_new = runge_kutta_4(u, t, dt, dx) # integration method, change manually---euler_f, runge_kutta_4
field[c, 0] = u_new
u = ti.cast(u_new, dtype=ti.f32)
t = ti.cast(t_new, dtype=ti.f32)
@ ti.func
def euler_f(u, t, dt, dx):
u_new = u + heat_equ(u, t, dx)*dt
t_new = t + dt
return u_new, t_new
@ ti.func
def runge_kutta_4(u, t, dt, dx):
k0 = heat_equ(u, t, dx)
k1 = heat_equ(u+(dt*k0/2.0), t+(dt/2.0), dx)
k2 = heat_equ(u+(dt*k1/2.0), t+(dt/2.0), dx)
k3 = heat_equ(u+dt*k2, t+dt, dx)
return u + dt*(k0+(2.0*k1)+(2.0*k2)+k3)/6.0, t+dt
@ ti.func
def d2udx2(u, dx):
d2u = ti.cast(dut, dtype=ti.f32)
for i in range(1, 31):
d2u[i] = (u[i+1]-2.*u[i]+u[i-1])/(dx**2.)
#print(d2u)
return d2u
@ ti.func
def heat_equ(u, t, dx, k=130.):
res = k*d2udx2(u, dx)
return res
field = ti.Vector.field(32, dtype=ti.f32, shape=(10000, 1)) # the compilation of all data points,
# in which the 'row' dimension must be
# the same as 'leng' parameter
ic = np.zeros(32)
ic[1:-1] = np.ones(30)
init_con = ti.Vector(list(ic), dt=ti.f32) # global initial condition
d = ti.types.vector(32, ti.f32) # global derivatives template
dut = d(np.zeros(32))
print(init_con)
inti_main()
#result = np.array(field)
print(field[9999, 0])