Home > K25104 > PDEs > stencil_order.m

stencil_order

PURPOSE ^

Ttaylor expands the numerical scheme prescribed by the stencil A.

SYNOPSIS ^

function [ final ] = stencil_order( A,order )

DESCRIPTION ^

 Ttaylor expands the numerical scheme prescribed by the stencil A.
 Input arguments:
   A, array holding the numerical stencil. The array should be a square 
   or cube of odd side length. The centre element will be interpreted
   as the centre of the stencil.
   order, highest power of h in the output
 Output arguments:
   final, string holding the expansion of the scheme

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ final ] = stencil_order( A,order )
0002 % Ttaylor expands the numerical scheme prescribed by the stencil A.
0003 % Input arguments:
0004 %   A, array holding the numerical stencil. The array should be a square
0005 %   or cube of odd side length. The centre element will be interpreted
0006 %   as the centre of the stencil.
0007 %   order, highest power of h in the output
0008 % Output arguments:
0009 %   final, string holding the expansion of the scheme
0010 
0011 [m,n,l]=size(A);  % find the size of A
0012 if m~=n;
0013     error('input array must be square or cube');
0014 elseif mod(n,2)~=1;
0015     error('matrix must have odd side length');
0016 elseif mod(order,1)~=0;
0017     error('input order must be an integer value');
0018 end
0019 
0020 c=(n+1)/2;      % find the centre index of the stencil
0021 syms Dx Dy Dz h;     % create symbolic variables Dx,Dy, Dz and h
0022 e=0;        % initilaize a symbolic variable to hold the error
0023 
0024 order = order+1;
0025 if l==1;    % if array is 2D
0026     for i=1:n;
0027         for j=1:n;
0028             tay=taylor(exp((i-c)*Dy*h+(j-c)*Dx*h),Dy,'order',order);
0029             tay=taylor(tay,Dx,'order',order);
0030             e=e+A(i,j)*tay;
0031         end
0032     end
0033 elseif l==n;    % if array is 3D
0034     for i=1:n;
0035         for j=1:n;
0036             for k=1:n;
0037             tay=taylor(exp((i-c)*Dy*h+(j-c)*Dx*h+(k-c)*Dz*h),Dy,'order',order);
0038             tay=taylor(tay,Dx,'order',order);
0039             tay=taylor(tay,Dz,'order',order);
0040             e=e+A(i,j,k)*tay;
0041             end
0042         end
0043     end
0044 else
0045     error('input array must be square or cube');
0046 end
0047 
0048 [cx,tx]=coeffs(e,h);
0049 cx=flip(cx);
0050 tx=flip(tx);
0051 
0052 final=blanks(1);
0053 hind=1;
0054 first=1;
0055 for k=1:order
0056     %temp=blanks(1);
0057     syms a;
0058     a=k;
0059     if hind<=length(cx);
0060         if isequaln(tx(hind),h^a);
0061             temp1=char(h^a);
0062             temp2=char(cx(hind));
0063             if first==1;
0064                 final=strcat(final,temp1,'{',temp2,'}');
0065                 first=0;
0066             else
0067                 final=strcat(final,' +',temp1,'{',temp2,'}');
0068             end
0069             hind=hind+1;
0070         end
0071     end
0072 
0073 end
0074 end
0075 
0076 
0077 
0078

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005