USER GUIDELINES OF: PARABOLIC PDE SOLVER Revision 2.11 of December 2, 2010 by M. MICCIO
MUC 1.0 ALLOWS TO SOLVE P.D.E. OF THE KIND : • SECOND ORDER; • PARABOLIC; • LINEAR; • COSTANT COEFFICIENTS.
INITIAL CONDITION : The initial condition can be a general function of x , not only a costant value u0 .
BOUNDARY CONDITIONS : The boundary conditions can be function of time: C(t) e F(t).
WHAT MUC 1.0 CAN DO: • Resolution of PDE with: explicit Euler method and implicit Crank&Nicholson method. • Check of the Mean Square Error between the solutions.
RESOLUTION OF P.D.E. (1) • Δ; • k . • Lenght of the domain (Lenght); • Number of points for discretization (# Points); • Time of last solution (Time); • Timestep for time discretization (Timestep); The parameters to indicate for each specific problem are :
RESOLUTION OF P.D.E. (2) It is necessary to specify initial and boundary conditions in the proper controls. Initial and boundary conditions can be function, respectevely, of x and t.
VALID FUNCTIONS FOR I.C. AND B.C. (1) abs(x) Absolute Value acos(x) Inverse Cosine acosh(x) Inverse Hyperbolic Cosine asin(x) Inverse Sine asinh(x) Inverse Hyperbolic Sine atan(x) Inverse Tangent atanh(x) Inverse Hyperbolic Tangent cos(x) Cosine cosh(x) Hyperbolic Cosine cot(x) Cotangent
VALID FUNCTIONS FOR I.C. AND B.C. (2) • csc(x) Cosecant (1/sin(x)) • exp(x) Exponential • expm1(x) Exponential (Arg - 1): (e^(x-1)) • getexp(x) Mantissa & Exponent (returns the exponent of x) • getman(x) Mantissa & Exponent (returns the mantissa of x) • int(x) Round To Nearest (rounds its argument to the nearest integer) • intrz(x) Round Toward 0 (rounds x to the nearest integer between x and zero) • ln(x) Natural Logarithm
VALID FUNCTIONS FOR I.C. AND B.C. (3) • lnp1(x) Natural Logarithm (Arg +1) • log(x) Logarithm Base 10 • log2(x) Logarithm Base 2 • max(x,y) Maximum • min(x,y) Minimum • mod(x,y) Quotient & Remainder • pow(x,y) x^y • rand( ) Random Number (0- 1) • rem(x,y) Remainder • sec(x) Secant [computes the secant of x, where x is in radians: (1/cos(x))] • sign(x) Sign (returns 1 if x is greater than 0, returns 0 if x is equal to 0, and returns -1 if x is less than 0)
VALID FUNCTIONS FOR I.C. AND B.C. (4) • sin(x) Sine • sinc(x) Sinc [computes the sine of x divided by x radians: (sin(x)/x)] • sinh(x) Hyperbolic Sin • sqrt(x) Square Root • tan(x) Tangent • tanh(x) Hyperbolic Tangent N.B.1 The independent variable can be indicated as either x or t N.B.2 The costant must be indicated as pi(1)
EULER METHOD (1) Timestep and spatial step must be chosen ( ) so that the Stability Parameter be: Otherwise the solution is unstable, and the warning led of the Stability Parameter lights on.
EULER METHOD (2) If the Stability Parameter is greater than 0.5 the solution is unstable. Pushing the button Euler Method Stability Safety, the s.p. is forced to 0.499 and the graph does not twist. Anyway the solution plotted is wrong if the Stability Parameter led is red.
CRANK-NICHOLSON METHOD • Implicit Method • Ever stable • Stability Parameter suggested:
Example of students’ test 1) Check the stability of the explicit method for the PDE having = 1and k = 0and further subject to a linear initialcondition: I.C.: u(x,0) = 0,4*x to a Dirichlet condition at the left boundary: B.C.1: u(0,t) =1 and a a mixedcondition at the right boundary: with: Lenght of the domain (Lenght) = 1 Number of points for discretization (# Points) = 40 Time of last solution (Time) = 0.3 Timestep for time discretization (Timestep) = 0,00033 2) How much is the spatial step? 3) What is the new value for the time step if the explicit method turns out unstable? 4) Discuss and comment the final diagram 5) Repeat integration with the Crank-Nicholson method and compare the results with the explicit method through the Mean Square Error
STABILITY EXAMPLE 1 – EULER METHOD Parameters: Lenght=1 #Points=40 Time=0,3 Timestep=0,00033 =1 k=0 I.C. u(x,o)=0,4x A=1 B=0 C=1 D=1 E=-1 F=0 Unstable Solution Initial Condition
STABILITY EXAMPLE 2 – EULER METHOD Parameters: Lenght=1 #Points=40 Time=0,02-0,03 Timestep=0,00055 =1 k=0 I.C. u(x,0)=sin(pi(1)*x) A=1 B=0 C=0 D=1 E=0 F=0
EXAMPLE 1 - CRANK & NICHOLSON METHOD Parameters: Lenght=1 #Points=40 Time=1,7 Timestep=0,00327024 =1 k=0 I.C. sin(pi(1)*x) A=1 B=0 C=abs(sin(pi(1)*x)) D=1 E=0 F=0
EXAMPLE 2 - CRANK & NICHOLSON METHOD Parameters: L=1 #Points=41 Time=1 Timestep=0,00082 =1 k=0 I.C. abs(sin(2*pi(1)*x)) A=1 B=0 C=abs(cos(4*pi(1)*x)) D=1 E=0 F=abs(sin(4*pi(1)*x))
COMPARING OF THE SOLUTIONS MUC 1.0 allows comparison of the solutions calculated with Euler and Crank-Nicholson methods through the function Mean Square Error.
COMPARING OF THE SOLUTIONS (2) The solution vectors –calculated by solving the same problem with each of the two methods – are compared with the MSE appearing in the proper indicator.
MEAN SQUARE ERROR Two vectors having n elements: can be compared calculating the MSE as:
MEAN SQUARE ERROR (2) The MSE found comparing two solutions can not indicate which solution is better: it only calculates a mean error between the solutions. To test the software, some numerical solutions have been compared with analytical solutions or commercial softwares. Here some of the results.
MEAN SQUARE ERROR (3) Parameters: Lenght=1; Time=0,05; =1; k=0; I.C. u(x,0)=0; A=1; B=0; C=1; D=1; E=0; F=1. 2 KINDS OF PROOF: Variable # Points and constant Timestep Constant # Points and variable Timestep
MUC 1.0 vs COMMERCIAL SOFTWARE It has been considered the problem of diffusion of a component through a layer of infinite lenght and unitary thickness.
FEM SOLUTION WITH A COMMERCIAL SOFTWARE AUTOMATIC MESH : Number of nodes = 16 Number of elements = 15 Max spatial step = 0,066667 Time = 0,1 Timestep = 0,001
MUC 1.0 SOLUTION Parametri: Lenght=1 # Points=16 Time=0,1 Timestep=0,01 =1 k=0 I.C. u(x,0)=0 A=1 B=0 C=1 D=1 E=0 F=1
MEAN SQUARE ERROR Comparing the Crank & Nicholson Solution Vector with the FEM solution for time 0,1 : Timestep=0,01 S.P.=1,1250 MSE=1,3 E-3 Timestep=0,001 S.P.=0,1125 MSE=3 E-5 Timestep=0,00001 S.P.=0,0013 MSE= 5,5 E-6
MUC 1.0 vs ANALYTICAL SOLUTION (1) Problem of monodimensional diffusion: Analytical solution:
MUC 1.0 vs ANALYTICAL SOLUTION (2) Lenght=1 # Points=16 Time=0,1 Timestep=0,0001
MUC 1.0 vs ANALYTICAL SOLUTION (3) MSEA/E = 5,3 E-7 MSEA/C&N = 7,3 E-7
Future developments: solving PDE with two and three spacial variables and systems of PDE. PARABOLIC PDE SOLVER by Ugo Avagliano student of Chemical and Food Engineering University of Salerno - Italy