Description of scalar

 SCALAR_HYP Numerical solution of scalar linear hyperbolic equations,
     formulas (3.7.1a), pag. 145, CHQZ2

      du/dt +beta du/dx=0           in (xa,xb) x (0,T). 
      initial cond.: u(x,t)=u0(x),  
      inflow cond.: u(xa,t)=ul(t)

  Space discretization: 1. Galerkin-Numerical Integration with LGL formulas,
                        2. Legendre collocation (strong form),
                        3. penalty,

  Time discretization:  RK4

  [x,u,err,Psi,Phi]=scalar_hyp(xa,xb,t0,T,beta,uex,u0,ul,nx,deltat,param);

 Input: xa,xb  = extrema of space domain Omega=(xa,xb)
        t0,T   = extrema of time domain [t0,T]
        beta   = transport coefficient
        uex    = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
        u0    = initial condition (u0=@(t)[u0(t)], with .*, .^, ./)
        ul    = inflow condition (ul=@(x)[ul(x)], with .*, .^, ./)
        nx     = spectral polynomial degree, the number of LGL points is
                 nx+1
        deltat = time step
        param(1)= space disc. scheme: 1=strong (3.7.2), pag. 146 CHQZ2, 
                                     2=GNI+ weak b.c. (3.7.5), pag. 147
                                     3=penalty (3.7.2a)+(3.7.9), pag. 148
        param(2)=tau : penalty parameter (used only if param(1)==3)

 Output: u = solution at final time T
         err= L_\infty norm error between numerical and exact solution at
                    final time T
         Psi = column array of evaluations Psi(t_n) (see def. 3.7.14),
               pag. 151, CHQZ2
         Phi = columen array of evaluations Phi(t_n) where
               Phi(t)=1/2 d/dt ||u||^2+1/2beta*(u^2(1)-u^2(-1))

 Possible input data:
 xa=-1;xb=1;
 beta=1.5;
 param(1)=2;param(2)=0.75
 nx=8;
 t0=0;T=4;
 deltat=1.e-4;

 Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
                    "Spectral Methods. Fundamentals in Single Domains"
                    Springer Verlag, Berlin Heidelberg New York, 2006.

0001 function [x,u,err,Psi,Phi]=scalar_hyp(xa,xb,t0,T,beta,uex,u0,ul,nx,deltat,param)
0002 % SCALAR_HYP Numerical solution of scalar linear hyperbolic equations,
0003 %     formulas (3.7.1a), pag. 145, CHQZ2
0004 %
0005 %      du/dt +beta du/dx=0           in (xa,xb) x (0,T).
0006 %      initial cond.: u(x,t)=u0(x),
0007 %      inflow cond.: u(xa,t)=ul(t)
0008 %
0009 %  Space discretization: 1. Galerkin-Numerical Integration with LGL formulas,
0010 %                        2. Legendre collocation (strong form),
0011 %                        3. penalty,
0012 %
0013 %  Time discretization:  RK4
0014 %
0015 %  [x,u,err,Psi,Phi]=scalar_hyp(xa,xb,t0,T,beta,uex,u0,ul,nx,deltat,param);
0016 %
0017 % Input: xa,xb  = extrema of space domain Omega=(xa,xb)
0018 %        t0,T   = extrema of time domain [t0,T]
0019 %        beta   = transport coefficient
0020 %        uex    = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
0021 %        u0    = initial condition (u0=@(t)[u0(t)], with .*, .^, ./)
0022 %        ul    = inflow condition (ul=@(x)[ul(x)], with .*, .^, ./)
0023 %        nx     = spectral polynomial degree, the number of LGL points is
0024 %                 nx+1
0025 %        deltat = time step
0026 %        param(1)= space disc. scheme: 1=strong (3.7.2), pag. 146 CHQZ2,
0027 %                                     2=GNI+ weak b.c. (3.7.5), pag. 147
0028 %                                     3=penalty (3.7.2a)+(3.7.9), pag. 148
0029 %        param(2)=tau : penalty parameter (used only if param(1)==3)
0030 %
0031 % Output: u = solution at final time T
0032 %         err= L_\infty norm error between numerical and exact solution at
0033 %                    final time T
0034 %         Psi = column array of evaluations Psi(t_n) (see def. 3.7.14),
0035 %               pag. 151, CHQZ2
0036 %         Phi = columen array of evaluations Phi(t_n) where
0037 %               Phi(t)=1/2 d/dt ||u||^2+1/2beta*(u^2(1)-u^2(-1))
0038 %
0039 % Possible input data:
0040 % xa=-1;xb=1;
0041 % beta=1.5;
0042 % param(1)=2;param(2)=0.75
0043 % nx=8;
0044 % t0=0;T=4;
0045 % deltat=1.e-4;
0046 %
0047 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0048 %                    "Spectral Methods. Fundamentals in Single Domains"
0049 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0050 
0051 %   Written by Paola Gervasio
0052 %   $Date: 2007/04/01$
0053 
0054 
0055 delt1=1/deltat;
0056 tt=(t0:deltat:T)';nt=length(tt);
0057 Phi=[]; UN1=[]; UL=[];Psi=[]; INTU=[];
0058 npdx=nx+1;
0059 [x,wx]=xwlgl(npdx); 
0060 [dx]=derlgl(x,npdx);
0061 jac=(xb-xa)*5.d-1;
0062 x=x*jac+(xa+xb)*5.d-1;
0063 alpha=25/12; alpha0=4; alpha1=-3; alpha2=4/3; alpha3=-0.25;
0064 
0065 
0066 tau=param(2)/(jac*wx(1));
0067 t=t0;
0068 un=zeros(npdx,1); un=u0(x);
0069 intu0=sum(un.*wx)*jac;
0070 UL=[UL;un(npdx)-ul(t)];
0071 INTU=[INTU;0];
0072 intu1ul=0;
0073 psi=0;
0074 Psi=[Psi;psi];
0075 
0076 
0077 
0078 
0079 % Construction of matrix A
0080 
0081 if param(1)==1
0082 % strong collocation (3.7.2)
0083 A=-beta*dx/jac;
0084 
0085 elseif param(1)==2
0086 % GNI+weak b.c. (3.7.5)
0087 A=beta*dx'*diag(wx);
0088 A(npdx,npdx)=A(npdx,npdx)-beta;
0089 
0090 elseif param(1)==3
0091 % penalty (3.7.2a) +(3.7.9)
0092 A=-beta*diag(wx)*dx;
0093 A(1,1)=A(1,1)-tau*beta*wx(1)*jac;
0094 
0095 end
0096 deltat2=deltat*0.5;
0097 
0098 
0099 
0100 
0101 %close all
0102 % fig=figure(...,
0103 %     'Name','Hyperbolic 1d problem',...
0104 %     'Visible','on')
0105 % axis([-1,1,-1,1]);
0106 % time loop explicit 4th order Runge Kutta
0107 for it=1:nt-1
0108 t=tt(it);
0109 if(param(1)==1)
0110 w1=A*un;
0111 w2=A*(un+deltat2*w1);
0112 w3=A*(un+deltat2*w2);
0113 w4=A*(un+deltat*w3);
0114 u=un+deltat/6*(w1+2*w2+2*w3+w4); u(1)=ul(tt(it+1));
0115 elseif(param(1)==2)
0116 w1=A*un; w1(1)=w1(1)+beta*ul(t); w1=(w1./wx)/jac;
0117 w2=A*(un+deltat2*w1); t=t+deltat2; w2(1)=w2(1)+beta*ul(t);w2=(w2./wx)/jac;
0118 w3=A*(un+deltat2*w2);w3(1)=w3(1)+beta*ul(t);w3=(w3./wx)/jac;
0119 w4=A*(un+deltat*w3); t=t+deltat2; w4(1)=w4(1)+beta*ul(t);w4=(w4./wx)/jac;
0120 u=un+deltat/6*(w1+2*w2+2*w3+w4);
0121 elseif(param(1)==3)
0122 w1=A*un; w1(1)=w1(1)+tau*beta*ul(t)*wx(1)*jac;w1=(w1./wx)/jac;
0123 w2=A*(un+deltat2*w1); t=t+deltat2; w2(1)=w2(1)+tau*beta*ul(t)*wx(1)*jac;
0124 w2=(w2./wx)/jac;
0125 w3=A*(un+deltat2*w2);w3(1)=w3(1)+tau*beta*ul(t)*wx(1)*jac;w3=(w3./wx)/jac;
0126 w4=A*(un+deltat*w3); t=t+deltat2; w4(1)=w4(1)+tau*beta*ul(t)*wx(1)*jac;
0127 w4=(w4./wx)/jac;
0128 u=un+deltat/6*(w1+2*w2+2*w3+w4);
0129 end
0130 
0131 intu=sum(u.*wx)*jac;
0132 t=tt(it+1); UL=[UL;u(npdx)-ul(t)];
0133 INTU=[INTU;intu-intu0];
0134 
0135 if(fix((it+1)/2)*2~=it+1)
0136 intu1ul=intu1ul+simpcz(tt(it-1:it+1),UL(it-1:it+1));
0137 psi=intu-intu0+beta*intu1ul;
0138 Psi=[Psi;psi];
0139 end
0140 nor=sum(u.*u.*wx)*jac;
0141 
0142 % plot(x,u);
0143 % title(['t=',num2str(t)]);
0144 % pause(0.01)
0145 %
0146 
0147 % evaluate Phi(t)=1/2 d/dt ||u||^2+1/2beta*(u^2(1)-u^2(-1))
0148 if(it==1)
0149 un=u; nor0=nor;
0150 elseif (it==2)
0151 un1=un; un=u;nor1=nor0; nor0=nor;
0152 elseif(it==3)
0153 un2=un1; un1=un; un=u;nor2=nor1;nor1=nor0; nor0=nor;
0154 elseif(it==4)
0155 un3=un2; un2=un1; un1=un; un=u;
0156 nor3=nor2;nor2=nor1;nor1=nor0; nor0=nor;
0157 else
0158 u1=u(1);u2=u(npdx);
0159 phi1=5.d-1*beta*(-u1^2+u2^2)+...
0160 5.d-1*delt1*(alpha*nor-(alpha0*nor0+alpha1*nor1+alpha2*nor2+alpha3*nor3));
0161 Phi=[Phi;phi1];
0162 un3=un2; un2=un1; un1=un; un=u;
0163 nor3=nor2;nor2=nor1;nor1=nor0; nor0=nor;
0164 end
0165 end
0166 
0167 
0168 % error  at  time T
0169 
0170 ue=uex(x,tt(nt));
0171 err=norm(u-ue);
0172 
0173 % fig=figure(...,
0174 %     'Name','Hyperbolic 1d problem Phi(t)',...
0175 %     'Visible','on')
0176 % plot(tt(5:nt-1),Phi);
0177 % ylabel('Phi(t)'); xlabel('t')

scalar_hyp

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE