求解杆的1-D热方程式在8次迭代后完全分解

Question

我正在使用有限的差异技术，或者实际上只是Euler的方法，但具有部分导数。在“水浴”中有两个端的标记边界条件。它在最初的7次迭代中衰减良好，但是在第八次，热信号无处可寻。它必须来自热WRT位置的第二部分衍生物，但我不知道原因。

import numpy as np
from math import *
import matplotlib.pyplot as plt

# Solving the Heat Equation d(u)/d(t) = k * d^2(u)/d(x)^2


tot=100 #just a number of points
L=10 #length of rod
tf=10 #final time
k=3 #k for konstant

x=np.linspace(0,L,tot) #1-D position
t=np.linspace(0,tf,tot) #time
u=np.zeros((tot,tot)) #heat(x,t)

dx=L/tot 
dt=tf/tot


for i in range(tot):
    u[i,0]=6*sin(pi/L*x[i]) #heat at t=0 boundary. Just some function of x



for i in range(tot): #heat at x=0 and x=L edge boundaries
    u[0,i]=0
    u[tot-1,i]=0

for i in range(tot): #i=time
    z=[] #z is just used for plotting
    z.append(u[0,i])
    for j in range(1,tot-1): #j=x-position
        u_xx=(u[j+1,i]-2*u[j,i]+u[j-1,i])/(dx**2) #get u_xx at j,i
        u_t=k*u_xx #solve eq for du/dt
        u[j,i+1]=u_t*dt + u[j,i] #integrate du/dt
        
        z.append(u[j,i])
    z.append(u[tot-1,i])
        
    plt.plot(z)
    plt.show()

    
#print(u[50])#iteration 50

Answer 1

数值方案可能是或不可能是可持续性。不稳定意味着解决方案中的任何小误差（例如，舍入错误）会很快增加。您使用的计划是不稳定的。对于线性方程，可以通过（例如）谐波分析进行方案稳定性分析。该分析表明，如果2 k dt/dx/dx＆lt; 1，欧拉方案将是稳定的。
无条件稳定方案（曲柄 - 尼古尔森计划）的例子如下：

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

def init(x):
    out=np.zeros(len(x))
    out[20:40]=2
    out[60:65]=4
    return out
def bound(t):
    return[np.zeros(t.shape),np.zeros(t.shape)]

def parab_eq(nt,nx,init,bound):
    def tridiag_alg(a,b,c,r):
        #a - first element not used
        #c - end element not used
        #b - main diagonal
        #a - lower diaganal
        #c - upper diaganal
        n=len(b)
        alf=np.zeros(n)
        beta=np.zeros(n)
        out=np.zeros(n)
        alf[0]=-c[0]/b[0]
        beta[0]=r[0]/b[0]
        for ii in range(1,n):
            alf[ii]=-c[ii]/(b[ii]+alf[ii-1]*a[ii])
            beta[ii]=(r[ii]-a[ii]*beta[ii-1])/(b[ii]+alf[ii-1]*a[ii])
        out[-1]=beta[-1]
        for ii in reversed(range(n-1)):
            out[ii]=alf[ii]*out[ii+1]+beta[ii]
        return out
    #nx,nt - number of intervals
    def solve_layer(u,coun_t):
        a=-kdif*dt/dx**2/2.0*np.ones(nx+1)
        b=(1+kdif*dt/dx**2.0)*np.ones(nx+1)
        c=-kdif*dt/dx**2/2.0*np.ones(nx+1)
        r=np.zeros(nx+1)
        for ii in range(1,nx):
            op1=kdif*dt/dx**2/2.0*(u[ii+1,coun_t]-2*u[ii,coun_t]+u[ii-1,coun_t])
            r[ii]=u[ii,coun_t]+op1
        r[1]+=-u[0,coun_t+1]*a[1]
        r[-2]+=-u[-1,coun_t+1]*c[-2]
        return tridiag_alg(a[1:-1],b[1:-1],c[1:-1],r[1:-1])

    t=np.linspace(0,tmax,nt+1)
    x=np.linspace(0,xmax,nx+1)

    dt=tmax/nt; dx=xmax/nx
    u=np.zeros((nx+1,nt+1))
    u[:,0]=init(x)
    u[0,:]=bound(t)[0]
    u[-1,:]=bound(t)[1]

    for count_t in range(nt):
        u[1:-1,count_t+1]=solve_layer(u,count_t)
    return t,x,u

Nx=100
K=2000
kdif=2
tmax=0.1
xmax=2

t,x ,result = parab_eq(K, Nx,init,bound)

f=plt.figure()
ax=plt.axes(xlim=(0,xmax),ylim=(-2,5))

tex = 'Numerical solution for: $\\frac{\\partial^2{u}}{\\partial{t^2}}=k\\frac{\\partial^2{u}}{\\partial{x^2}}

plt.title(tex, fontsize=20, color='black')

ax.grid()
plt.ylabel('u')
plt.xlabel('x')
line,=plt.plot([],[],'k')
def init():
    line.set_data([],[])
    return line
def update(frame):
    line.set_data(x,result[:,frame])
    return line

ani = FuncAnimation(f,update, init_func=init,frames=range(0,K,2), interval=100,repeat=True)
plt.show()