STAG_SCALAR_HYP Numerical solution of scalar linear hyperbolic equations,
formulas (3.7.1a), pag. 145, CHQZ2 by staggered grids method (pag. 149,
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: Staggered grids method (pag. 149, CHQZ2)
Time discretization: RK4
[x,u,err,Psi,Phi]=stag_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
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;
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]=stag_scalar_hyp(xa,xb,t0,T,beta,uex,u0,ul,nx,deltat) 0002 % STAG_SCALAR_HYP Numerical solution of scalar linear hyperbolic equations, 0003 % formulas (3.7.1a), pag. 145, CHQZ2 by staggered grids method (pag. 149, 0004 % CHQZ2) 0005 % 0006 % du/dt +beta du/dx=0 in (xa,xb) x (0,T). 0007 % initial cond.: u(x,t)=u0(x), 0008 % inflow cond.: u(xa,t)=ul(t) 0009 % 0010 % Space discretization: Staggered grids method (pag. 149, CHQZ2) 0011 % 0012 % Time discretization: RK4 0013 % 0014 % [x,u,err,Psi,Phi]=stag_scalar_hyp(xa,xb,t0,T,beta,uex,u0,ul,nx,deltat,param); 0015 % 0016 % Input: xa,xb = extrema of space domain Omega=(xa,xb) 0017 % t0,T = extrema of time domain [t0,T] 0018 % beta = transport coefficient 0019 % uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./) 0020 % u0 = initial condition (u0=@(t)[u0(t)], with .*, .^, ./) 0021 % ul = inflow condition (ul=@(x)[ul(x)], with .*, .^, ./) 0022 % nx = spectral polynomial degree, the number of LGL points is 0023 % nx+1 0024 % deltat = time step 0025 % 0026 % Output: u = solution at final time T 0027 % err= L_\infty norm error between numerical and exact solution at 0028 % final time T 0029 % Psi = column array of evaluations Psi(t_n) (see def. 3.7.14), 0030 % pag. 151, CHQZ2 0031 % Phi = columen array of evaluations Phi(t_n) where 0032 % Phi(t)=1/2 d/dt ||u||^2+1/2beta*(u^2(1)-u^2(-1)) 0033 % 0034 % Possible input data: 0035 % xa=-1;xb=1; 0036 % beta=1.5; 0037 % nx=8; 0038 % t0=0;T=4; 0039 % deltat=1.e-4; 0040 % 0041 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0042 % "Spectral Methods. Fundamentals in Single Domains" 0043 % Springer Verlag, Berlin Heidelberg New York, 2006. 0044 0045 % Written by Paola Gervasio 0046 % $Date: 2007/04/01$ 0047 0048 0049 delt1=1/deltat; 0050 tt=(t0:deltat:T)';nt=length(tt); 0051 Phi=[]; UN1=[]; UL=[];Psi=[]; INTU=[]; 0052 npdx=nx+1; 0053 [x,wx]=xwlgl(npdx); 0054 [dx]=derlgl(x,npdx); 0055 [xg,wxg]=xwlg(nx); 0056 jac=(xb-xa)*5.d-1; 0057 alpha=25/12; alpha0=4; alpha1=-3; alpha2=4/3; alpha3=-0.25; 0058 0059 % eta = interpolation matrix from LG grid (xg) to LGL grid (x) 0060 0061 [eta]=intlag_lg(xg, wxg, x); 0062 0063 % ltrm = interpolation matrix from LGL grid to LG grid 0064 0065 [ltrm]=intlag_lgl(x, xg); 0066 0067 % dg Legendre Gauss derivative matrix 0068 0069 [dg] = derlg(xg,nx); dg=dg/jac; 0070 0071 0072 x=x*jac+(xa+xb)*5.d-1; 0073 xg=xg*jac+(xa+xb)*5.d-1; 0074 0075 0076 t=t0; 0077 un=zeros(nx,1); un=u0(xg); 0078 intu0=sum(un.*wxg)*jac;u2=eta(npdx,:)*un; 0079 UL=[UL;u2-ul(t)];INTU=[INTU;0]; 0080 intu1ul=0; 0081 psi=0; 0082 Psi=[Psi;psi]; 0083 0084 0085 utilde=zeros(npdx,1); 0086 0087 0088 deltat2=deltat*0.5; 0089 0090 for it=1:nt-1 0091 t=tt(it); 0092 % 0093 utilde=eta*un; 0094 utilde(1)=ul(t); 0095 % 0096 w1=-beta*dg*(ltrm*utilde); 0097 0098 % 0099 utemp=un+deltat2*w1; 0100 utilde=eta*utemp; 0101 t=t+deltat2; utilde(1)=ul(t); 0102 w2=-beta*dg*(ltrm*utilde); 0103 0104 % 0105 utemp=un+deltat2*w2; 0106 utilde=eta*utemp; 0107 utilde(1)=ul(t); 0108 w3=-beta*dg*(ltrm*utilde); 0109 0110 % 0111 utemp=un+deltat*w3; 0112 utilde=eta*utemp; 0113 t=t+deltat2; utilde(1)=ul(t); 0114 w4=-beta*dg*(ltrm*utilde); 0115 0116 u=un+deltat/6*(w1+2*w2+2*w3+w4); 0117 nor=sum(u.*u.*wxg)*jac; 0118 t=tt(it+1); u1=ul(t);u2=eta(npdx,:)*u; 0119 intu=sum(u.*wxg)*jac; 0120 UL=[UL;u2-u1];INTU=[INTU;intu-intu0]; 0121 if(fix((it+1)/2)*2~=it+1) 0122 intu1ul=intu1ul+simpcz(tt(it-1:it+1),UL(it-1:it+1)); 0123 psi=intu-intu0+beta*intu1ul; 0124 Psi=[Psi;psi]; 0125 end 0126 0127 0128 % evaluate Phi(t)=1/2 d/dt ||u||^2+1/2beta*(u^2(1)-u^2(-1)) 0129 if(it==1) 0130 un=u; nor0=nor; 0131 elseif (it==2) 0132 un1=un; un=u;nor1=nor0; nor0=nor; 0133 elseif(it==3) 0134 un2=un1; un1=un; un=u;nor2=nor1;nor1=nor0; nor0=nor; 0135 elseif(it==4) 0136 un3=un2; un2=un1; un1=un; un=u; 0137 nor3=nor2;nor2=nor1;nor1=nor0; nor0=nor; 0138 else 0139 phi1=5.d-1*beta*(-u1^2+u2^2)+... 0140 5.d-1*delt1*(alpha*nor-(alpha0*nor0+alpha1*nor1+alpha2*nor2+alpha3*nor3)); 0141 Phi=[Phi;phi1]; 0142 un3=un2; un2=un1; un1=un; un=u; 0143 nor3=nor2;nor2=nor1;nor1=nor0; nor0=nor; 0144 end 0145 end 0146 0147 % compute error at T 0148 ue=uex(xg,tt(nt)); 0149 err=norm(u-ue); 0150 0151 return