Heat diffusion

From https://en.wikipedia.org/wiki/Heat_equation

"In mathematics and physics (more specifically thermodynamics), the heat equation is a parabolic partial differential equation. The theory of the heat equation was first developed by Joseph Fourier in 1822 for the purpose of modeling how a quantity such as heat diffuses through a given region. Since then, the heat equation and its variants have been found to be fundamental in many parts of both pure and applied mathematics". 

The idea is to see how Rebol or Red can solve this problem. Let’s imagine a coin heated at its center at T0.0 and calculate the temperature distribution toward the edges of the coin. We get a nice Gaussian curve.

Now, let’s look at how the phenomenon changes over time (from 0.0 to 0.10). We still see nice Gaussian curves whose mean decreases over time. We can also model the cooling of the metal part. 

Due to vector datatype  (https://github.com/Oldes/Rebol3) these equations are very easy to solve.

#!/usr/local/bin/r3

REBOL [

]

b2d: import 'blend2d ;--use blend2d (draw module)

cv2: import 'opencv ;--use OpenCV for visualisation

;--Global Physical Parameters

k: 0.5 ;--Diffusion coefficient (Thermal diffusivity coefficient of material)

l: 1.0     ;--Domain size

itime: 0.1   ;--Integration time

;--Global Numerical Parameters

nx: 100     ;--Number of grid points

nt: 1000   ;--Number of time steps

dx: l / (nx - 1)   ;--Grid step (Spatial step size: 0.0101010101010101)

dt: itime / nt   ;--Grid step (Time step size :0.0001)

;--Based on Oldes's code

linSpace: func [

    "Generates n linearly spaced numbers from a to b (inclusive)."

    a [number!]  "Start"

    b [number!]  "End"

    n [integer!] "Number of samples"

    /no-end      "Don't include end value in the result"

    /local vec div step i ;--local values

][

    vec: make vector! [f64! :n]

    div: either no-end [n][n - 1]

    step: (b - a) / div

    repeat i n [vec/:i: a + (step * (i - 1))]

    vec

]

;--Initialisation (using global variables)

init: does [

x: linSpace 0.0 1.0 nx ;--linear spacing between 0.0 and 1.0

t: (x * 2.0 * pi) ;--t is a vector

repeat i nx [t/:i: sin t/:i] ;--in radians

;--use blend2d commands in plot block

plot: compose [

font %/System/Library/Fonts/Geneva.ttf

fill black

text 0x10  10 "+1.0"

text 0x204 10 "0"

text 0x390 10 "-1.0"

pen black

line 5x200 410x200 

]

]

;--Heat diffusion equation (using global variables)

diffusion: does [

vect: make vector! [f64! :nx]

posy: 0

n: 0

while [n <= nt] [

j: 2

while [j < nx] [

vect/:j: dt * k * (t/(j - 1) - (2.0 * t/:j) + t/(j + 1)) / (dx ** 2)

++ j

]

repeat j nx [t/:j: t/:j + vect/:j]

;--Plot every 100 time steps

if zero? n % 100 [ 

acolor: random white

posy: posy + 15

either n = 0 [tt: 0.0] [tt: round/to n % (n * dt) 0.001]

tt: ajoin ["time: " tt]  

pos: as-pair 340 posy

append plot reduce ['fill (acolor) 'text (pos) 12 (tt) 'pen (acolor) 'line]

xx: 2

foreach val t [

append plot as-pair xx * 4 200 - (val * 195) 

++ xx 

]

]

++ n

]

]

;************************ Main Program ************************

init 

diffusion

img: b2d/draw 410x400 :plot

img: cv2/resize img 150% 

cv2/imshow/name :img "Heat Diffusion"

print "Any key to close"

cv2/waitKey 0