FTCS-Algorithmus für die Wärmegleichung

Question

Ich versuche, den FTCS-Algorithmus für die 1-dimensionale Wärmeleitungsgleichung in Python zu implementieren.

import numpy as np 

L = 1 #Length of rod in x direction 
k = 0.3 #Thermal conductivity of rod 
tmax = 5 #how many seconds 
nx = 100 #number of spacial steps 
nt = 100; #number of time steps 
xi = np.linspace(0,L,nx) 
ti = np.linspace(0,tmax,nt) 
dx = L/(nx-1) 
dt = tmax/(nt-1) 
r = k*dt/(dx)**2 
r2 = 1-2*r 

u=np.zeros((nt,nx)) 
#IC 
phi = 100; 
for x in range(0,nx):  
    u[0][x] = phi 

#BC 
for t in range(0,nt):  
    u[t][0] = 0; 
    u[t][nx-1] = 0 
#FTCS Algorithm 

for t in range(0,nt-1): #timestep 
    for x in range(1,nx-2): 
     u[t+1][x] = r*(u[t][x-1]+ u[t][x+1]) + r2*(u[t][x]) 

aber ich bin nicht plausible Werte für u bekommen [t] [x] = u (x, t) unter Berücksichtigung meiner Anfangsbedingungen von 100. dh sie sprengen und geben Sie mir dumme Werte wie ‚4.11052068e + 221 'gibt es eine schlechte Programmierpraxis, an der ich teilnehme, die den Algorithmus zerstört? Oder habe ich den Algorithmus gerade falsch implementiert?

EDIT: Ich denke, ich habe herausgefunden, dass es ist, weil der Algorithmus stabil ist, wenn und nur wenn r < 1/2. Die Zahlen explodieren einfach, weil mein r etwa 2,5 oder so ist, aber wenn jemand andere Fehler sehen kann, lass es mich wissen !!

Answer 1

Original Question

Die Stabilität der FTCS Regelung auf der Größe der konstanten r Scharniere. Wenn , dann werden Rundungsfehler, die bei jedem Schritt eingeführt werden, exponentiell abnehmen. Wenn r>1/2, dann werden diese Rundungsfehler exponentiell zunehmen. (Wie Sie in Ihrer Bearbeitung angedeutet haben).

Klein-ish Fehler

1. dx = L/nx und dt = tmax/nt. (Wenn das für Sie verwirrend ist, stellen Sie sich den Fall von nx = 2 und L = 1 ... dann dx = 0.5 vor).
2. L und sollte Floats sein. (dh - L = 1.0).
3. Sie aktualisieren die Temperatur u nicht an der vorletzten Position auf der Stange. Stellen Sie sicher, dass sich die for loop over-Position auf nx-1 und NICHT nx-2 erstreckt.

Things

betrachten

4. Sie könnten mit realen Werten/Einheiten für Zeit, Länge, Wärmeleitfähigkeit und Ausgangstemperatur berücksichtigen.

5. Wenn Sie die Randbedingungen (unter dem Kommentar #BC) gelten, wird nur die Zeit t=0 hat Werte von u, die nicht Null sind, so müssen Sie alle Zeitschritte nicht durchlaufen.

6. Während es notwendig ist, dass Sie durch die Zeit iterieren (daher die Zeit for Schleife), können Sie tatsächlich die räumliche Ableitung vektorisieren. Sie können Ihre for-Schleife durch folgende ersetzen: for t in range(0, nt-1): u[t+1, 1:nx-1] = r*(u[t,0:nx-2] + u[t,2:]) + r2*(u[t,1:nx-1])