<jupyter_start><jupyter_text>One-dimensional diffusion equationv1.71, 14 February 2018, by <NAME> Demonstrates the numerical solution of the one-dimensional diffusion equation.The time required for the forecast to run is limited by the deliberate delay in the animation.  So the effectiveness using numpy arrays, rather than lists, is not obvious yet. In future notebooks, numpy arrays will be essential for the forecast to run quickly.  So take special note of the numpy forecast.<jupyter_code>#import matplotlib
#matplotlib.rcParams.update({'font.size': 22})
from IPython.core.display import HTML
from IPython.display import display,clear_output
import time as Time
import math
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import urllib.request
HTML(urllib.request.urlopen('http://metrprof.xyz/metr4323.css').read().decode())
#HTML( open('metr4323.css').read() ) #or use this, if you have downloaded metr4233.css to your computer<jupyter_output><empty_output><jupyter_text>Check out the wikipedia: [Heat equation in one dimension](http://en.wikipedia.org/wiki/Heat_equationDerivation_in_one_dimension) Let $u(x,t)$ be *temperature* and $q(x,t)$ be *heat flux density*.  You would probably be happier with $T$ as the symbol for temperature, but that's the way it goes...$$q = -k \frac{ \partial u}{ \partial x}$$$$\rho c_p \frac{ \partial u}{ \partial t} = - \frac{ \partial q}{ \partial x}$$If we let $$F \equiv \frac{q}{\rho c_p}$$then $F$ will play the role of a *temperature flux density*. Our coupled PDEs now appearas:$$F = -\frac{k}{\rho c_p} \frac{ \partial u}{ \partial x}$$$$\frac{ \partial u}{ \partial t} = - \frac{ \partial F}{ \partial x}$$And we can further tidy up the notation with:$$\alpha  \equiv \frac{k}{\rho c_p}$$We can combine the coupled first order equations to obtain the familiar one-dimensional diffusion equation:$$\frac{ \partial u}{ \partial t} = \alpha \frac{\partial^2 u} {\partial x^2}$$ Re-run the model from here<jupyter_code># make the grid
N = 101 # number of x grid points for u
xmax = 20 # 0 <= x <= xmax
dx = xmax/(N-1.) # grid width
xu = [dx*i for i in range(N)] # grid for u
xF = [.5*dx + dx*i for i in range(N-1)] # grid for F, between u
# initialize u, F and alpha
init = 'gauss'
#init = 'square'
#init = 'sines'
#init = 'twodeltax'
xc = xmax/2. # center of domain
if init == 'gauss':
    ui = [ math.exp( -(x-xc)**2 ) for x in xu ] # Gaussian for initial field of u
elif init == 'square':
    ui = [1 if abs(x-xc)<2.001 else 0 for x in xu] # square wave
elif init == 'sines':
    tp = 2*math.pi/xmax
    ui = [.5+.25*math.sin(tp*x)+.25*math.sin(3*tp*x) for x in xu]
else:
    ui = [.5+.25*math.cos(tp*x)+.25*math.cos(51*tp*x) for x in xu] 
#    print("not implemented")

Fi = [0]*(N-1) # flux is set to zero, will be computed    

# a quick look at ui
plt.plot(xu,ui)
plt.ylim((0,1.1))
# set physical parameter:
alpha=1. # diffusion coefficient

# set some crunching parameters

# choose one value for crunch:

crunch='py'  # python lists
#crunch='lc' # python list comprehensions
#crunch = 'np' # python numpy arrays

dt = .01 # try dt=.01, .0201, .022, .03
tstop = 4. # stop when t>tstop
dplot = .1 # time increment between plots
# run the model forward in time
# initialize fields:
if crunch == 'np':
    u = np.array(ui) # a numpy array
    F = np.array(Fi) 
    dudt = np.zeros(N)
else:
    u = ui[:] # this slice copies the entire list
    F = Fi[:] 
    dudt = [0.]*N  # list of N zeros
    
t = 0. # initial value for t    
tplot = 0. # time t for next plot

f, ax = plt.subplots()
display(f)

while t < tstop + dt/2.:
    
    if tplot-dt/2.<t<tplot+dt/2.: #plot
        Time.sleep(1./10.) #Slows down the computation
        ax.clear()
        ax.plot(xu,u)
        ax.text(1,.9,crunch,fontsize=22)
        ax.text(1,.5,'t={0:5.3f}'.format(t),fontsize=22)
        clear_output(wait=True)
        ax.set_xlabel('x')
        ax.set_ylabel('u',rotation=0)
        ax.set_ylim((-.2,1.2))
        display(f)
        tplot = t+dplot
        if t < dt/2.: Time.sleep(2.)
        if t > tstop: break
            
    # step fields forward in time with Euler method.  Nothing fancy!
    if crunch=='py': # use python lists
        for n in range(N-1):
            F[n] = alpha*(u[n] - u[n+1])/dx
        for n in range(N-2):
            dudt[n+1] = (F[n] - F[n+1])/dx
        for i in range(N):
            u[i] = u[i] +dt*dudt[i] 
            
    elif crunch=='lc': # use python lists, with  list comprehensions
        F = [ alpha*(u[n] - u[n+1])/dx for n in range(N-1) ]
        dudt = [0.]+ [ (F[n] - F[n+1])/dx for n in range(N-2) ] + [0.]
        u = [ u[i] + dt* dudt[i] for i in range(N) ]
        
    elif crunch=='np': # use numpy arrays
        F = alpha*( u[:-1] - u[1:] )/dx         
        dudt[1:-1] = ( F[:-1] - F[1:] )/dx
        u = u + dudt*dt 
        
    t = t + dt
    
plt.close() 
beta = 4*alpha*dt/(dx*dx)
print('beta was:', beta )<jupyter_output>beta was: 0.9999999999999998
{'fraction_non_alphanumeric': 0.15287263961430295, 'fraction_numerical': 0.027721976697468863, 'mean_word_length': 2.932859399684044, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '249363', 'license': 'MIT', 'n_tokens_mistral': 1889, 'n_tokens_neox': 1779, 'n_words': 595, 'path': 'Fiedler_revised/n040_diffusionpde1d.ipynb', 'repo_name': 'ATMOcanes/Fiedler_notebooks', 'size': '4985'}