### Numerical Method

 The diffusion simulation is base on the rules : 1. The temperature changes follow the formula of diffusion            Ut = D * Uxx 2. The boundary condition is isolated            Uxn+1 = UxnHere is the core part of this program.   The formula of diffusion will be convert into program code.   The Formula of Diffusion : Ut = D * UxxU is a function that record the temperature values alone x-axis.    Meanwhile , itchanges by time t.That means :        U = U(x,t)Ux means the slope of U at x point.   Or ,how temperature changes when it moves alone x-axis at x point.Let's consider about temperature value function U at point x , U(x). Assume h is very small distance alone x-axis , or h -> 0        2h  =  U(x+1)  -  U(x)  =  U(x)  -  U(x-1)        U(x+2h) = U(x+1)Then by definition        Ux(x) = [ U(x+h) - U(x-h) ] / 2h        Uxx(x) = [ Ux(x+h) - Ux(x-h) ] / 2h        Ux(x+h) = [ U(x+h+h) - U(x+h-h) ] / 2h = [ U(x+1) - U(x) ] / 2h         Ux(x-h) = [ U(x-h+h) - U(x-h-h) ] / 2h = [ U(x) - U(x-1) ] / 2h Then we have ,           Uxx(x) = [ Ux(x+h) - Ux(x-h) ] / 2h = [ U(x+1) - 2U(x) + U(x-1) ] / (4h^2)Apply this to y-axis         Uyy(y) = [ U(y+1) - 2U(y) + U(y-1) ] / (4h^2) Apply this to Ut (For 2D , dt -> 0 ),        Ut = D * ( Uxx + Uyy) =  D * { [ U(x+1) - 2U(x) + U(x-1) ] / (4h^2) + [ U(y+1) - 2U(y) + U(y-1) ] / (4h^2) }             = D * [ U(x+1,y) + U(x-1,y) + U(x,y+1) + U(x,y-1) - 4U(x,y) ] / (2h^2)        U(t+1) = U(t) + dt * D1 / h^2 * [ U(x+1,y) + U(x-1,y) + U(x,y+1) + U(x,y-1) - 4U(x,y) ]This is numerical form of Diffusion Equation ( Ut = D * Uxx ) !Each calculation will get the temperature value after time interval dt.This method is called Finite Difference Method.[The Meaning of dt * D1 / h^2 ]Assume c = dt * D1 / h^2                   = dt * D2 / h1^2                    = dt * Dcof / dx^2 Here h1 = 2h = U(x+1) - U(x) = dx = The width of each calculation nod        D2 = Dcof = The diffusion coefficient of a materialThe unit of :        U(x) : Celcius        dt : sec        dx : cm        Dcof : cm^2 / secThe dt , dx should be very small.    Then we can set         dt = 0.01 sec        dx = 0.01 cmSo , we can use 500 nods to present a material that has 5cm width.And each calculation will get the result of the time passed by 0.01 sec In this program , the c = 0.200 , it is a constant in here.        c = dt * Dcof / dx^2        Dcof = 0.002 (cm^2/sec)This says : IF , we set dt = 0.01 sec , dx = 0.01 cmThen , the program shows the simulation result of a material that has Dcof = 0.002 (cm^2/sec)[Remark]You'll get wrong result if writing code like this :        U(x) = U(x) + c * [ U(x-1) - 2U(x) + U(x+1) ]This is the correct form :         U2(x) = U(x) + c * [ U(x-1) - 2U(x) + U(x+1) ]And , the value of c should not exceed 0.5.[Hint About Programming]class Material {  bool bActive;   float value;} u[n];In this 1D case , the user can pick up any point from 1 to (n-1) for diffusion.Say , n == 100 , but the material only occupy u[5] to u[10]So , // initializefor(int i = 5 ; i<=10 ; i++ ) {   u[i].bActive = true;  u[i].value = ini_value ;}// for calculationfor(int i = 1 ; i<=(n-1) ; i++) {  if( u[i].bActive == true ) {     // deal with boundary condition     if( u[i-1].bActive == false ) u[i-1].value = u[i].value;      if( u[i+1].bActive == false ) u[i+1].value = u[i].value;     normal calculation   }}// draw on screenshow( u[i] );The example on the web page only calculate 50 times per sec , update show( ) once after more calculations can make it much faster. Copyright © 2010 bdragon All rights reserved.bdragong@gmail.comWebRepOverall rating