Solving the Heat Diffusion Equation (1D PDE) in Python

[@BelmanCinematography]

import numpy as np import matplotlib.pyplot as plt L = 0.1 #Grosor de la pared en metros n = 10 #Numero de nodos utilizados T0 = 20 #Temperatura inicial T1s = 40 #Temperatura superficial en cara 1 T2s = 20 #Temperatura superficial en cara 2 dx = L/n alpha = 0.000000054713 #Difusividad termica K/(Rho*C_p) t_final = 1800 #Tiempo final en segundos dt = 60 x = np.linspace(dx/2, L-dx/2, n) T = np.ones(n)*T0 dTdt = np.empty(n) t = np.arange(0, t_final, dt) for j in range(1,len(t)): plt.clf() for i in range(1, n-1): dTdt[i] = alpha*(-(T[i]-T[i-1])/dx**2+(T[i+1]-T[i])/dx**2) dTdt[0] = alpha*(-(T[0]-T1s)/dx**2+(T[1]-T[0])/dx**2) dTdt[n-1] = alpha*(-(T[n-1]-T[n-2])/dx**2+(T2s-T[n-1])/dx**2) T = T + dTdt*dt plt.figure(1) plt.plot(x,T) plt.axis([0, L, 0, 50]) plt.xlabel('Distance (m)') plt.ylabel('Temperature (C)') plt.show() plt.pause(0.01)

[@vsl8902]

I'm looking for a numerical code of Fourier equation. KZbin algorithm bring me here. This is amazing! Thanks

[@lauraburgstaller8682]

do you know, how to save this animation as .mp4?

[@riyazhudda3772]

Excellent explanation. Thank you! I would replace all the dT/dt calculations as follows: # define this once discrete_calc = lambda Ti_1,Ti,Ti1 : alpha * (-(Ti-Ti_1)/dx**2 + (Ti1-Ti)/dx**2) Now replace the definition of dT/dt vector everywhere dTdt[i] = discrete_calc(T[i-1],T[i],T[i+1]) dtDt[0] = discrete_calc(T1s, T[0], T[1]) dtDt[n-1] = discrete_calc(T[n-2],T[n-1],T2s)

[@sidebo1]

This is a great video, with one small "insidious" bug (meaning it's a bug that introduces an error, but won't cause program to crash, and in this case, would not introduce much error). It's with the 2 boundary conditions, where the distance from the node to the surface should be dx/2, not dx. It's an easy fix in those 2 lines of code.

[@PRACHISINGH-en4yz]

Can you please let me know at what line, this change needs to be done?

[@manuelo1738]

Hello Mr. Powell, I would like to ask you a question. One the boundary conditions shound't the distance be dx/2? Also, on the first for loop the end should be len(t)+1 to finish in the last time step. Kind regards

[@boudzzzzsssdddfffg]

For those struggeling with getting the animation right, it might help to add the command plt.ion() in the first for loop. That should stop it from opening the tab every run.

[@sayanjitb]

Yes, it just activates an interactive environment.

[@PRACHISINGH-en4yz]

Thank you so much. Your comment saved me!

[@JustMoseyinAround]

Very well done sir. Thank you

[@stefanvet8135]

Nice and clear explanation

[@bryanvazquez6422]

thank's a lot for this video!, i've got a question... how did you make to "animate" the profile temperature at 23:47?

[@shivamm9882]

That's what he told in the video. The reason to put plt.pause was to reduce the speed of the graph generation.

[@alip6001]

plt.show()

[@gabry88chesterfield]

plt.ion() before the for cycle!

[@ShashwatSharan]

@@gabry88chesterfield finally it's moving! thanks!

[@biopellet]

really helpful stuff

[@KV-gy2mr]

This is so cool. Could you also show us how to do this by solving a system of equations using linalg i.e implicit finite difference method?

[@aurelia8028]

Yeah I'd like to see how that is done aswell. I couldn't quite wrap my head around the wikipedia article with matrices and all that jazz and how I'd code that

[@TriThom50]

Is there a way to extend this to be able to apply heat at different points along the x direction? Also, would this result still hold if this was a thin rod, meaning could we consider a thin rod using the 1D equation too? Because in one of the lectures I know you mentioned a wall.

[@camimartinez6963]

thank you i love you you saved my life

[@JamesVestal-dz5qm]

The heat transfer equation is from advanced transport phenomena 2, and were coding in python in reaction engineering. Brandon Tatum messaged me what does chemical engineering have to do with my presidential campaign and honestly i dont know. Brandon, my goal is to get healthy as a chemical engineer. Im learning about chemical engineering while waiting for my presidency to happen.

[@Lonkines]

Thank you very much for the video!

[@m35926]

Quick question, what numerical method technique should I use if I want to find heat flow through a series of flat plates? So imagine we take your picture and keep the heat flow the same direction but push the wall on it's side so its layers are stacked up like pancakes. Any help with that? And awesome video.

[@sarazahoor3110]

I tried to run your code and it is giving me the following error, $python main.py File "main.py", line 25 dTdt [0] = alpha*(-(T[0]-T1s)/dx**2+(T[1]-T[0])/dx**2) ^ SyntaxError: invalid syntax

[@sayanjitb]

Can you post the whole program that you have written?

[@ecarlson]

Fantastic tutorial!

[@LateMax2359]

I wrote my own (very similar) program independently before discovering your guide. I've found that both our programs have the same bug. This may be a conceptual error with how I understand the heat transfer equation, but as dx decreases, I would expect the simulation to get more accurate, and not do anything weird. But, as I lower the thickness of the material or increase the number of nodes I find that the numbers begin blasting off by several orders of magnitudes immediately. I'm not sure why this occurs, but I'm (marginally) happier knowing that this isn't a unique problem. My current method is by lowering the number of nodes as much as is reasonable, and if continuing to lower it would be unreasonable I assume that steady state heat transfer begins almost immediately. I'd be interested to know what your opinion on the issue is.

[@bluesybluesko9922]

I might be late to the party, but here it goes : you have a so called CFL condition that connects to both dx and dt, and that condition should be lower than one, otherwise you start getting "junk" results and the method becomes unstable

[@LateMax2359]

@@bluesybluesko9922 yes, very late, though very much appreciated

[@gluijk]

Very didactic example. It can be easily replicated in R using vector notation (hyper fast update since the 'for i' spatial loop becomes unnecesary): T=matrix(T0, nrow=1, ncol=n) T[1]=Tleft T[n]=Tright indices=which(col(T)!=1 & col(T)!=n) x=seq(from=dx/2, to=L-dx/2, length.out=n) # position array val=alpha*(dt/dx^2) for (j in 0:N) { plot(x, T, type='l', col='red', lwd=2, xlim=c(0,L), ylim=c(min(Tleft,Tright), max(Tleft,Tright)), ylab='T (ºC)', main=paste0('Iteration: ', j, '/',N, ', t=', j*dt, 's')) abline(v=c(x[1], x[n]), lty=2) abline(h=0) # Iterate T in vector notation T[indices]=T[indices]+val*(T[indices+1] - 2*T[indices] + T[indices-1]) }

[@nabeelkubba]

I'm trying to run this, but the figure only shows one profile then stops, after I close it, it then shows the next profile in time, it is not flowing as a video! Any ideas what is going on? I am using Python 3.7 and attempting to run it from Visual Studio. Thanks!

[@luchtbuks1]

Use plt.draw() instead of plt.show()

[@AnupumPant]

Amazing! Works like a charm. I added a gaussian source in the middle too. But can you explain how can I add a zero flux boundary condition.

[@euyin77]

dTdx = 0 at x = L

[@AdityaRoyiseast]

hey, why did you use linspace in 14 and arange in 19? why does error pops when i interchange them?

[@BelmanCinematography]

Instead of an animation it gives me a list of images at the end, I'm using Python via Jupyter 5.5.0, does anyone know how to turn it into an animation?

[@ariaalizadeh8067]

I had the same problem. I just had to bring "plt.show()" outside of the for loop!

[@juandadamo]

jupyter.brynmawr.edu/services/public/dblank/jupyter.cs/Examples/Animations%20in%20Matplotlib.ipynb A working example.

[@dshshajskj6143]

CAn you make one using fem?

[@mariospapanicolaou4631]

THANKS A LOT FRIEND!!!!

[@alip6001]

sorry can you explain more about T=T+dTdt*dt, because my phthon says, operands could not be broadcast together with shapes (600,) (10,)

[@andriesesun3073]

excuse me,I meet the same question. Could you tell how to solve this question?

[@laraibquamar7156]

sir could you please explain ADI scheme for 2D heat equation in python.

[@joaopedronevesgoldenstein8566]

Did you get it Laraib? I need it too