NONPERIODIC_BURGERS Numerical solution of non periodic Burgers equation
\partial u/\partial t +u \partial u/\partial x -\nu \partial^2 u/\partial x^2
=0 in Omega, \forall t>t_0
u(x,t_0)=u_0(x) in Omega
+ b.c.
Exact solution is given by formula (3.1.8), pag. 119, CHQZ2
Space discretization: Galerkin-Numerical Integration with LGL
formulas.
Time discretization: Explicit Euler, RK2, RK4
[u,err]=nonperiodic_burgers(xa,xb,t0,T,visc,nx,deltat,tscheme,bc)
Input: xa,xb = extrema of space domain Omega=(xa,xb)
t0,T = extrema of time domain [t0,T]
visc = viscosity
nx = spectral polynomial degree, the number of LGL points is
nx+1
deltat = time step
tscheme= time discretization scheme: 1=EE, 2=RK2, 3=RK4
bc = choice of boundary conditions: 1 == Dirichlet
2 == Neumann
Output: u = solution at final time T
err= L_\infty norm error between numerical and exact solution at
final time T
Possible input data:
xa=-1;xb=1;
visc=0.01;
tscheme=4;bc=1;
nx=16;
t0=0;T=1;
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 [u,err]=nonperiodic_burgers(xa,xb,t0,T,visc,nx,deltat,tscheme,bc) 0002 % NONPERIODIC_BURGERS Numerical solution of non periodic Burgers equation 0003 % 0004 % \partial u/\partial t +u \partial u/\partial x -\nu \partial^2 u/\partial x^2 0005 % =0 in Omega, \forall t>t_0 0006 % u(x,t_0)=u_0(x) in Omega 0007 % + b.c. 0008 % 0009 % Exact solution is given by formula (3.1.8), pag. 119, CHQZ2 0010 % Space discretization: Galerkin-Numerical Integration with LGL 0011 % formulas. 0012 % Time discretization: Explicit Euler, RK2, RK4 0013 % 0014 % [u,err]=nonperiodic_burgers(xa,xb,t0,T,visc,nx,deltat,tscheme,bc) 0015 % 0016 % Input: xa,xb = extrema of space domain Omega=(xa,xb) 0017 % t0,T = extrema of time domain [t0,T] 0018 % visc = viscosity 0019 % nx = spectral polynomial degree, the number of LGL points is 0020 % nx+1 0021 % deltat = time step 0022 % tscheme= time discretization scheme: 1=EE, 2=RK2, 3=RK4 0023 % bc = choice of boundary conditions: 1 == Dirichlet 0024 % 2 == Neumann 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 % 0030 % 0031 % Possible input data: 0032 % xa=-1;xb=1; 0033 % visc=0.01; 0034 % tscheme=4;bc=1; 0035 % nx=16; 0036 % t0=0;T=1; 0037 % deltat=1.e-4; 0038 % 0039 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0040 % "Spectral Methods. Fundamentals in Single Domains" 0041 % Springer Verlag, Berlin Heidelberg New York, 2006. 0042 0043 % Written by Paola Gervasio 0044 % $Date: 2007/04/01$ 0045 0046 0047 % Definition of exact solution and first derivative to set Neumann b.c 0048 0049 Uex=['1+(x-t)/(t+1).*(4/sqrt(t+1)*exp(-(x-t).^2/(0.04*(t+1))))',... 0050 './(1+4/sqrt(t+1)*exp(-(x-t).^2/(0.04*(t+1))))']; 0051 Uex1=['-4*(50*(t+1)^(1/2)*t^2-100*x*t*(t+1)^(1/2)-4*exp(-25*(x-t)^2/(t+1))*t',... 0052 '+50*x^2*(t+1)^(1/2)-(t+1)^(1/2)*t-(t+1)^(1/2)-4*exp(-25*(x-t)^2/(t+1)))',... 0053 '*exp(-25*(x-t)^2/(t+1))/((t+1)^(1/2)+4*exp(-25*(x-t)^2/(t+1)))^2/(t+1)^2']; 0054 uex=vectorize(Uex); 0055 uex1=vectorize(Uex1); 0056 uex=@(x,t)[eval(uex)]; 0057 uex1=@(x,t)[eval(uex1)]; 0058 0059 0060 0061 % space discretization 0062 0063 npdx=nx+1; 0064 [x,w]=xwlgl(npdx); 0065 dx=derlgl(x,npdx); 0066 jac=(xb-xa)*5.d-1; 0067 x=x*jac+(xb+xa)*5.d-1; 0068 A=stiff_1d_sp(w,dx,jac); 0069 A=visc*A; 0070 0071 t=t0; 0072 u0=uex(x,t); 0073 0074 fig=figure(... 0075 'Name','Non periodic Burgers solution',... 0076 'Visible','on'); 0077 title(['t=',num2str(t)]); 0078 %U1=uex(x); 0079 %plot(x,u0,x,U1,'r',x,zeros(npdx,1),'+'); 0080 %pause(0.001) 0081 k=0; 0082 0083 % time loop 0084 0085 while t< T 0086 if tscheme==1 0087 % Explicit Euler in time 0088 u=ee(t,deltat,@fburgers,u0,dx,A,w,visc,uex,uex1,bc); 0089 elseif tscheme==2 0090 % RK2 in time 0091 u=rk2(t,deltat,@fburgers,u0,dx,A,w,visc,uex,uex1,bc); 0092 elseif tscheme==4 0093 % RK4 in time 0094 u=rk4(t,deltat,@fburgers,u0,dx,A,w,visc,uex,uex1,bc); 0095 end 0096 t=t+deltat; 0097 k=k+1; 0098 U1=uex(x,t); 0099 % if rem(k,100)==0 0100 % err=norm(u(1:nx)-U1(1:nx),inf); 0101 % fprintf('t=%13.6e, nx=%d, err_inf=%13.6e, tscheme=%d \n ',t,nx, err,tscheme) 0102 % end 0103 plot(x,u,x,U1,'r',x,ones(npdx,1),'+'); 0104 axis([-1,1,0.8,1.2]) 0105 title(['t=',num2str(t)]); 0106 pause(0.001) 0107 u0=u; 0108 if (sum(isnan(u))>1) 0109 fprintf('Nan solution: instability due to too large time-step') 0110 err=100; 0111 return 0112 end 0113 % end time loop 0114 end 0115 U1=uex(x,t); 0116 err=norm(u(1:nx)-U1(1:nx),inf); 0117 fprintf('t=%13.6e, nx=%d, err_inf=%13.6e, tscheme=%d \n ',t,nx, err,tscheme) 0118 0119 0120 %U1=uex(x,t); 0121 %err=norm(u(1:nx)-U1(1:nx),inf); 0122 %disp('N, error') 0123 %[nx, err] 0124 plot(x,u,x,U1,'r',x,ones(npdx,1),'+'); 0125 axis([-1,1,0.8,1.2]) 0126 title(['t=',num2str(t)]); 0127 %pause(0.001) 0128 0129 return