rkdg_o3.m

用有限元分析的法解决一类偏微分方程

%RKDG with order 3 in both time and space,for solving Ut+Ux=0,with U(x,0)=sin(x),basis bases on the legendre polynomials
for k=1:4
    Nx=20*(2^(k-1))
    dx=2*pi/Nx;dt=0.05/(2^(k-1));T=6;Nt=T/dt;x=zeros(Nx+1,1);U0=zeros(3,Nx+1);U1=U0;CFL=dt/dx;
  % initial value
  for j=2:Nx+1
      x(j)=(j-3/2)*dx;
      U0(1,j)=(cos(x(j)-dx/2)-cos(x(j)+dx/2))/dx;
      U0(2,j)=3*(-cos(x(j)-dx/2)-cos(x(j)+dx/2)+2*(sin(x(j)+dx/2)-sin(x(j)-dx/2))/dx)/dx;
      U0(3,j)=(45/4)*((8/(dx^2)-2/3)*(cos(x(j)+dx/2)-cos(x(j)-dx/2))+4*(sin(x(j)+dx/2)+sin(x(j)-dx/2))/dx)/dx;
  end
      U0(:,1)=U0(:,Nx+1);
  % matrix
  MA=1/dx*[-1,-1,-2/3;3,-3,-2;-15/2,15/2,-5];MA2=1/dx*[1,1,2/3;-3,-3,-2;15/2,15/2,5];
  %RKDG with order 3
   for n=1:Nt
     for j=2:Nx+1
         Utemp1(:,j)=U0(:,j)+dt*(MA*U0(:,j)+MA2*U0(:,j-1));
     end
         Utemp1(:,1)=U0(:,1)+dt*(MA*U0(:,1)+MA2*U0(:,Nx));
         Utemp2(:,1)=3*U0(:,1)/4+Utemp1(:,1)/4+dt*(MA*Utemp1(:,1)+MA2*Utemp1(:,Nx))/4;
     for j=2:Nx+1
         Utemp2(:,j)=3*U0(:,j)/4+Utemp1(:,j)/4+dt*(MA*Utemp1(:,j)+MA2*Utemp1(:,j-1))/4;
         U1(:,j)=U0(:,j)/3+2*Utemp2(:,j)/3+2*dt*(MA*Utemp2(:,j)+MA2*Utemp2(:,j-1))/3;
     end
      U1(:,1)=U1(:,Nx+1);U0=U1;
  end
  % L2 error 
  syms y;
  errorsquare=0;
  for j=2:Nx+1
      errorsquare=errorsquare+int((U0(1,j)*1+U0(2,j)*2*(y-x(j))/dx+U0(3,j)*(4*((y-x(j))/dx)^2-1/3)-sin(y-T))^2,y,x(j)-dx/2,x(j)+dx/2);
  end
  error(k)=subs(sqrt(errorsquare))
end %end the cycle of k
for k=2:4
    order(k)=log(error(k-1)/error(k))/log(2);
end
error=subs(error)
order=subs(order)