"""
heat.py
"""
#! /usr/bin/python
# -*- coding: utf-8 -*-
# interest of this program :
# numerical experiment of heat equation (1D)
# especially time step (dt) stability
doc = """
# Heat equation
# physics : d phi / d t = alpha nabla^2 phi
# alpha = 1
# numeric : (phi_1 - phi_0) / dt = nabla2(phi_0)
# phi_1 = phi_0 + dt * nabla2(phi_0)
"""
import sys
import os
import math
import cmath
import matplotlib.pyplot
import numpy
x_min = -float(10)
x_max = float(10)
x = [ ]
nb_phi0_plot = 0
phi_initial = [ ]
def step(n):
"""
# état initial :
# definition on xmin < 0 < xmax
# xmin < x < xmin/2 : phi(x) = 0
# xmin/2 < x < xmax/2 : phi(x) = 1 # îlot central à 1
# xmax/2 < x < xmax : phi(x) = 0
"""
global x, phi_initial
if nb_phi0_plot == 0 :
print("step function phi_0(x) on [0.5 * x_min ; 0.5 * x_max ]")
x = [ float(0) for i in range(n) ]
phi_0 = [ float(0) for i in range(n) ]
for i in range(n) :
x[i] = x_min + (x_max - x_min) * float(i) / float(n-1)
if 0.5 * x_min < x[i] and x[i] < 0.5 * x_max :
# phi = step
phi_0[i] = float(1)
phi_initial = list(phi_0)
return (phi_0)
def window(i, phi) :
# we add a x symmetry on the 2 boundaries
# as if the system was surrounded by same systems
n = len(x)
x0 = x[i]
if i < n-1 :
x1 = x[i+1]
else : # i == n-1
x1 = 2 * x[i] - x[i-1]
dx = x1 - x0 # constant : does not need to be centered
# x0 boundary : phi(x0-a) = phi(x0+a)
iminus1 = i-1
if iminus1 < 0 : iminus1 +=2
iplus1 = i+1
if iplus1 > n-1 : iplus1 -=2
phim1 = phi[iminus1]
phi0 = phi[i]
phip1 = phi[iplus1]
return (phim1, phi0, phip1, dx)
def nabla2(phi) :
""" d^2phi / dx^2 : (1D laplacian)
# dphi[i+1/2] = (phi[i+1] - phi[i]) / dx
# dphi[i-1/2] = (phi[i] - phi[i-1]) / dx
# d2phi[i] = (dphi[i+1/2] - dphi[i-1/2]) / dx
# d2phi[i] = (phi[i+1] - 2 phi[i] + phi[i-1]) / dx^2
"""
n = len(x)
d2_phi = [ float(0) for i in range(n) ]
for i in range(n) :
# with x symmetry boundary conditions :
(phim1, phi0, phip1, dx) = window(i, phi)
d2_phi[i] = (phip1 - float(2) * phi0 + phim1) / pow(dx, 2)
return d2_phi
def draw_phi_initial( ):
matplotlib.pyplot.plot(x, phi_initial)
def plot_phi(x, phi) :
# displays the position of phi versus x
# print("phi=", phi)
matplotlib.pyplot.text( 0. , 0. , "O")
matplotlib.pyplot.plot(x, phi)
def mean(phi) :
# boundary points belong to two domains
phi_moy = float(0)
n = len(phi)
phi_moy += 0.5 * abs(phi[0])
for i in range(1, n-1) :
phi_moy += abs(phi[i])
phi_moy += 0.5 * abs(phi[n-1])
phi_moy /= (n-1)
return phi_moy
def analysis(phi) :
phi_moy = mean(phi)
phi_initial_moy = mean(phi_initial)
if phi_moy > float(2) * phi_initial_moy :
print(" phi_moy=", phi_moy, "(ref :", phi_initial_moy,")")
print(" -> Warning : function phi diverges !")
def heat(n, dt) :
"""
# étude de la convergence de la solution du problème de la chaleur
# selon la discrétisation dx en espace
# selon la discrétisation dt en temps
# notion : pas de temps limite de stabilité dt = 0.5 dx²
"""
global nb_phi0_plot
# print(doc)
dt_stability = 0.5 * (x_max -x_min)**2 / float((n-1)**2)
if dt > dt_stability :
print("# Pas de temps trop grand (divergence) :")
elif dt == dt_stability :
print("# Pas de temps limite de la stabilité (oscillations) :")
elif dt < dt_stability :
print("# Pas de temps assurant la stabilité de l'algorithme :")
print(" n =", n, "; dt =", dt, " (dt_stability =", dt_stability, ")")
temps = float(0)
# initial state
step(n)
# initial state drawn once
if nb_phi0_plot == 0 :
matplotlib.pyplot.title("t = 0")
plot_phi(x, phi_initial)
nb_phi0_plot += 1
# heat evolution :
phi_0 = list(phi_initial)
for i in range(1000) :
d2_phi = nabla2(phi_0)
# print("nabla^2 phi : d2_phi=", d2_phi)
phi_1 = [ float(0) for j in range(n) ]
for j in range(n) :
# phi_1 = phi_0 + dt * nabla2(phi_0)
phi_1[j] = phi_0[j] + dt * d2_phi[j]
phi_0 = phi_1
# print("phi_0[",point,"]=", phi_0[point])
if i < 10 or \
i < 100 and i % 20 == 19 :
# matplotlib.pyplot.title("t = " + str(i+1))
plot_phi(x, phi_0)
# final state
# detection of divergence : negative phi
analysis(phi_0)
title = "t = " + str(i+1) + "(n=" + str(n) + ", dt=" + str(dt) + ")"
matplotlib.pyplot.title(title)
plot_phi(x, phi_0)
draw_phi_initial()
matplotlib.pyplot.show( )
# heat equation
# equation de la chaleur
# relation de stabilite
# stability relation : dt < 0.5 dx^2
help(heat)
for n in [ 11, 21, 31 ] :
dx = (x_max -x_min) / float(n-1)
dt_stability = 0.5 * dx**2
heat(n, 1.1 * dt_stability) # diverges
heat(n, 1.0 * dt_stability) # limit diverges (does not converge)
heat(n, 0.9 * dt_stability) # stable mais oscillant
heat(n, 0.5 * dt_stability) # stable
"""
"""