%------------------------------------------------------------------% % SHO_eigen_C1_L2.m % % Finite Element Eigensolution of 1-D % % Time Independent Schrodinger Equation (TISE) for SHO % % % % By Ed Akin, Rice University, as modification of Example 8.10.2 % % of Y.W. Kwon "The Finite Element Method Using Matlab", CRC, 1996 % % % % Problem description: - Y,xx + (Pot - 2E) Y = 0, with % % cubic C1 Hermite polynomials, with Y and Y,x at each of 2 nodes. % % % % Here: Pot is from the Simple Harmonic Operator (SHO), Pot = x^2 % % % % Variable descriptions % % p = element potential matrix % % k = element stiffness matrix, kk = system stiffness matrix % % m = element mass matrix, mm = system mass matrix % % index = equation numbers associated with each element % % bcdof = equation numbers associated with boundary conditions % % x_e = element nodal coordinates % %------------------------------------------------------------------% format long % need high precision % USE CONTROL nel = 10; % number of elements (then try 20) show = 4; % number of eigenvalues to show list = 0; % number of eigenvectors to list plots = 2; % number of waveforms, expected values to plot odd = 0; %1; % flag odd or even waveforms tleng = 5; % total (half)region length, approximates inf % END CONTROL nnel = 2; % number of nodes per element ndof = 2; % number of dofs per node (for C1) nnode = (nnel-1)*nel+1; % total number of nodes in system sdof = nnode*ndof; % total system dofs, = 22 here leng = tleng/nel; % size of TISE elements % specify Dirichlet bc equation numbers (1 fake eigenvalue each) if ( odd ) % Odd solutions: Y(0)=0, Y(x_max) = 0 bcdof (1) = 1; % Y at node 1 is constrained bcdof (2) = sdof-ndof+1; % Y at last node is constrained fprintf('SHO for odd modes with %g cubic elements \n', nel) else % Even solution: Y,x(0)=0, Y(x_max) = 0 bcdof (1) = 2; % Y,x at node 1 is constrained bcdof (2) = sdof-ndof+1; % Y at last node is constrained fprintf('SHO for even modes with %g cubic elements \n', nel) end % if bc type kk = zeros (sdof,sdof); % initialize system stiffness matrix mm = zeros (sdof,sdof); % initialization of system mass matrix index = zeros (nel*ndof,1); % initialization of index vector % define the Gaussian quadrature data on -1 to +1 pt = [ -0.906179845938663992797627 -0.538469310105683091036314 ... 0.000000000000000000000000 0.538469310105683091036314 ... 0.906179845938663992797627 ]; wt = [ 0.23692688505618908751426 0.47862867049936646804129 ... 0.56888888888888888888889 0.47862867049936646804129 ... 0.23692688505618908751426 ]; for iel = 1:nel % loop for the total number of elements index = feeldof1 (iel,nnel,ndof); % get element equation numbers x_e (1) = leng*(iel-1); % x left x_e (2) = leng*iel; % x right [k,m] = myTISE_C1_L2 (x_e,pt,wt); % compute C1 element matrices kk = feasmbl1 (kk,k,index); % assemble element into system mm = feasmbl1 (mm,m,index); % assemble element into system end % for iel elements [kn,mn] = feaplycs (kk,mm,bcdof); % apply the boundary conditions [Vecs, Diag] = eig (kn,mn); % get eigenvectors, eigenvalues for j = 1:sdof % extract diagonal matrix E (j) = Diag (j,j); % (fake bc values = 12345) end % for j [E, Order] = sort (E); % reorder eigenvalues, lo to hi % list the results if ( show > 0 ) for j = 1:show % list eigenvalues only fprintf ('Eigenvalue %g = %g \n', j, E(j)) % values end % for j end % if show if ( list > 0 ) for j = 1:list % list eigen value-vector pairs fprintf ('Eigenvalue %g = %g \n', j, E(j)) % values fprintf ('Eigenvector:') % vectors disp(Vecs(:,Order(j))) end % for j end % if list %_____________________________________________________________________ % Odd exact results are 2E = 3, 7, 11, 15, 19, ... % Here nel=10, tleng=5 give 3.00002, 7.00026, 11.0014, 15.0092, ... % Even exact results are 2E = 1, 5, 9, 13, 17 % Here nel=10, tleng=5 give 1, 5.00009, 9.00062, 13.0033, ... % Increasing nel and/or tleng increases accuracy (& problem size) % ================================================================== % post-process for ,,, and related graphs % by rebuilding element arrays since kn, mn destroyed by eig, and % = Y'*mn*Y, = Y'*kn*Y also equal the element sums of % Y_e'*m*Y_e and Y_e'*k*Y_e, etc if ( plots > 0 ) % then post-process x = [0:(nnode-1)]*leng; % x-axis for plots Y_e = zeros (nnel*ndof,1); % pre-allocate element dof wave = zeros (sdof,1); % pre-allocate waveform dof YY = zeros (nel,1); YHY = zeros (nel,1) ; YxY = zeros (nel,1) ; for p = 1:plots % consider each mode clf % clear plot frame hold on; grid ; % extract this waveform, get expected values from elements wave = Vecs(:, Order(p)); % waveform for this mode for iel = 1:nel % loop for the total number of elements x_e (1) = leng*(iel-1); % x left x_e (2) = leng*iel; % x right [k,m] = myTISE_C1_L2 (x_e,pt,wt); % re-compute C1 matrices [X] = my_exp_C1_L2 (x_e,pt,wt); % expectation matrices index = feeldof1 (iel,nnel,ndof); % element dofs numbers Y_e = wave (index); % el dof via vector subscripts YY(iel) = Y_e'*m*Y_e; % should sum to unity YHY(iel) = Y_e'*k*Y_e; % expected Hamiltonian YxY(iel) = Y_e'*X*Y_e; % expected x end % for iel fprintf('Plot eigenvalue %g = %g \n', p, E(p)) fprintf('Element contributions to : \n') disp(YY) Y_Y = sum(YY); Y_H_Y = sum(YHY); Y_x_Y = sum(YxY); fprintf(' = %g \n', Y_Y); % normalization fprintf(' = %g \n', Y_H_Y); % domain integral fprintf(' = %g \n', Y_x_Y); % domain integral fprintf('Expected 2E = %g \n', Y_H_Y/Y_Y); % expected E %b fprintf('Expected x = %g \n', Y_x_Y/Y_Y); % SHO expected x == 0 % scale function value, not slopes too ymax = max(max([wave(1:2:sdof),wave(1:2:sdof).^2])); axis([0,tleng,min(wave(1:2:sdof)),ymax]) if ( odd ) title(['SHO Piecewise cubic waveform for odd mode ', ... int2str(p), ', 2E = ', num2str(E(p))]) else title(['SHO Piecewise cubic waveform for even mode ', ... int2str(p), ', 2E = ', num2str(E(p))]) end % if even-odd xlabel(['X-axis (atomic units), for ', int2str(nel), ' elements']) ylabel('Node Value - o (solid), Y ^2 (dashed)') %b plot (x,wave(1:2:sdof),'ko'); %b plot (x,wave(1:2:sdof),'k-'); % waveform %b plot (x,wave(1:2:sdof).^2,'k--'); % Y^2 opt = 3; % plot choice plot_cubic (nel,x,wave,opt); % opt 1-Y 2-Y,x 3-Y & Y^2 % print hardcopy, pause, go to next mode print ('-dpsc', ['waveform_', int2str(p), '_graph']) hold off; v_text = ['Created waveform_', int2str(p), '_graph.ps']; fprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) pause (4) ; % note: could plot piecewise constant YY here for each el over x end % for p plots end % if post-proces %_____________________________________________________________________ % Related references: % 1. J.E. Akin, "Finite Element Analysis with Error Estimates" % www.owlnet.rice,edu/~mech517 (under Books link) % 2. L.R. Ram-Mohan, "Finite Element and Boundary Element Applications % in Quantum Mechanics", Oxford Univ. Press, 2002 % 3. W. Schweizer, "Numerical Quantum Dynamics", Kluwer, 2001 % end SHO_eigen_C1_L2