rkldg_o3.m

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

%RKLDG with order 3 in both time and space,for solving Ut=Uxx,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.0016/(2^(2*k-2));T=2;Nt=T/dt;x=zeros(Nx+2,1);U0=zeros(3,Nx+2);U1=U0;V=U0;CFL=dt/(dx^2);
  % 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);U0(:,Nx+2)=U0(:,2);
% matrix
  MA1=1/(dx^2)*[9,-9,6;21,-21,14;22.5,-22.5,15];MA2=1/(dx^2)*[9,7,2;-27,-21,-6;67.5,52.5,15];
  MA=1/(dx^2)*[-18,2,-8;6,-54,16;-90,60,-90];
  %RKLDG with order 3
   for n=1:Nt
     for j=2:Nx+1
         Utemp1(:,j)=U0(:,j)+dt*(MA1*U0(:,j+1)+MA*U0(:,j)+MA2*U0(:,j-1));
     end
         Utemp1(:,Nx+2)=Utemp1(:,2);Utemp1(:,1)=Utemp1(:,Nx+1);
     for j=2:Nx+1
         Utemp2(:,j)=3*U0(:,j)/4+Utemp1(:,j)/4+dt*(MA1*Utemp1(:,j+1)+MA*Utemp1(:,j)+MA2*Utemp1(:,j-1))/4;
     end
         Utemp2(:,Nx+2)=Utemp2(:,2);Utemp2(:,1)=Utemp2(:,Nx+1);
     for j=2:Nx+1
         U1(:,j)=U0(:,j)/3+2*Utemp2(:,j)/3+2*dt*(MA1*Utemp2(:,j+1)+MA*Utemp2(:,j)+MA2*Utemp2(:,j-1))/3;
     end
         U1(:,Nx+2)=U1(:,2);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)-exp(-2)*sin(y))^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)