source: trunk/sample/miscellaneous/FDM/FDM.mso @ 977

#*---------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2016 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2016 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*----------------------------------------------------------------------
* Author: Argimiro R. Secchi
* COPPE/UFRJ - Group of Modeling, Simulation, Control,
*                               and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#

using "types";

Model MDF
        
        PARAMETERS
        opcao as Switcher (Valid=["MDF1"],Default="MDF1");
        N as Integer (Brief="Number of discretization points");
        xi as Real(Brief="Lower boundary");
        xf as Real(Brief="Uppper boundary");

        VARIABLES
        x(N+1) as Real(Brief="mesh points");
        dx as Real(Brief="mesh interval");

        y(N+1) as Real(Brief="Dependent variable");
        dif1x(N+1) as Real(Brief="First derivative of y");
        dif2x(N) as Real(Brief="Second derivative of y"); 
        
        EQUATIONS
        
        #mesh uniform interval
        dx=(xf-xi)/N; 

        #mesh points
        for i in [1:N+1] do
                
                x(i) = xi + (i-1)*dx;
        
        end

        switch opcao
        
                case "MDF1":
                        #First derivative approximation at boundary x=x0
                        dif1x(1)=(y(2)-y(1))/(dx);
        
                        #First derivative approximation at boundary x=xf
                        dif1x(N+1)=(y(N+1)-y(N))/(dx);
        
                        #First derivative approximation at internal points 
                        for i in [2:N] do
                                dif1x(i)=(y(i+1)-y(i-1))/(2*(dx));
                        end
        
                        dif2x(1) = 0; # dummy
                        #Second derivative approximation at internal points 
                        for i in [2:N] do
                                dif2x(i)=(y(i+1) - 2*y(i) + y(i-1))/((dx)^2);
                        end
                
                #case "MDF2":
                
                        #ADD OTHER SCHEMES
                        
        end             

end