heat equation with nonflux boundary condition

16 Views Asked by At

I am trying to solve this following PDE

u_t=u_{xx}

with nonflux boundary condition such that $u_x=0$. I create artificial solution to test my solver where the true solution should be

u(x,t)=cos(\pi x/L)\exp{-(\pi/L)^2t}

This is my solver where I ensured the CFL condition. I used explicit method and centered for spatial points. But when I plot the error, it does not follow the 1st order of convergence. I tried both vectorized and loop.

clear all;close all;
    L=5;
    Time=.5;
    
    
    CFL=0.5;
    
    h=[1,0.5,0.1,0.01];
    error = zeros(1,length(h));
    
    for i=1:length(h)
    %     dx=1;
        dx=h(i);
        dt=dx^2*CFL;
        
        k=(pi/L)^2;
        
        heat = @(xx,tt) cos(pi.*xx./L).*exp(-k.*tt);
        initial = @(xx) cos(pi.*xx./L);
        
        x=0:dx:L;
        t=0:dt:Time;
        N=length(t);
        M=length(x);
        
        [X,T] = meshgrid(x,t);
        mu=dt/dx^2;
        
        data=heat(X,T);
        u0 = initial(x);
        
        
        solu=zeros(N,M);
        solu(1,:)=u0;
        
        for tstep=1:N-1
            un=zeros(1,M)';
            uc = solu(tstep,:)';
        
    %         un(2:end-1) = uc(2:end-1) + mu * ...
    %             (uc(1:end-2) - 2.0 * uc(2:end-1) + uc(3:end));
    % 
    %         un(1)=un(2);
    %         un(end)=un(end-1);
    
            for xstep=2:M-1
                un(xstep)=uc(xstep)+mu*( uc(xstep-1) - 2*uc(xstep) + uc(xstep+1) );
            end
            un(1)=un(2);
            un(end)=un(end-1);
            solu(tstep+1,:) = un;
            plot(un);
            hold on
        end
        error(i)=norm(data-solu,inf);
    end
    
    % figure; % Opens a new figure window
    % surf(X, T, sol-data);
    % 
    % % Adding labels and title
    % xlabel('space');
    % ylabel('time');
    % zlabel('Z-axis');
    % title('Surface Plot of Data Matrix');
    %%
    loglog(h,h.^2*error(1))
    hold on;
    loglog(h,h*error(1)*1.7)
    hold on;
    loglog(h,error)
    legend('2nd','1st - 1.7','error')
0

There are 0 best solutions below