contour function to show temperature on a 2d plate at 0, 1, 2, 5, 10, 100 seconds for discretised pde with jacobi solver. implicit euler for time and central differencing for space. using pointers in matlab

To plot the contour of the temperature on the 2D plate at specific times with the discretised PDE using Jacobi's solver, Implicit Euler for time and central differencing for space in MATLAB, we assume a rectangular plate over the domain (xmin, xmax)= (0,1) and (ymin, ymax)= (0,1) with N discretisation points for both x and y. We have to solve the heat diffusion equation in 2-D in the form:

The initial condition of the above heat diffusion equation is given by:

The boundary conditions for the above equation are:

Here, we use Jacobi Iterations to solve the diffusion equation at each time step. Given the Jacobi's method to numerically solve this problem, we only pass pointers as Jacobi's function input. A function jacobi is used to compute Jacobi's method with implicit Euler for time, central differences for space.

After getting the temperature matrix (T) at different times, we create the contour plots using the matlab function contourf.

Here's the MATLAB code for the above simulation.

main.m
%Define constants 
N = 50; %number of discretisation points
tmax = 100; %maximum time
dt = 1; %time step
dx=1/N;
dy=1/N;
k=1;
alpha=k*dt/(dx^2);

%Define initial temperature, T at t=0
T_old=zeros(N+1);
T_new=zeros(N+1);
for i=1:N+1
    for j=1:N+1
        if (i>1 && i<N+1 && j>1 && j<N+1)
            T_old(i,j)=sin(pi*(i-1)*dx)*sin(pi*(j-1)*dy);
        else
            T_old(i,j)=0;
        end
    end
end

%Create the grids 
[X,Y]=meshgrid(0:dx:1,0:dy:1);

%Create the figure
figure();

%Loop through time 
for t=0:dt:tmax
    
    %Copy the T_old to new array
    T_new(:,:)=T_old(:,:);

    %Jacobi Iterations
    for iter=1:500
        
        %Iterate over space (Exclude boundary points)
        for i=2:N
            for j=2:N
                T_old(i,j)=0.25*(T_new(i+1,j)+T_new(i-1,j)+T_new(i,j+1)+T_new(i,j-1))-alpha*T_new(i,j);
            end
        end
    end
    
    %Update T at each time step
    T_old = T_new;

    %Plot the contourf
    if t==0 || t==1 || t==2 || t==5 || t==10 || t==100
        contourf(X,Y,T_new,20,'LineStyle','none')
        title(['Temperature Contour at t = ',num2str(t)])
        colorbar
        pause(0.1)
    end
end
1171 chars
57 lines

Here, we first initialise the constants N, tmax, dt, dx, dy and k with discretisation and thermal properties of the plate.

Next, we define the initial temperature on the 2D plate at t=0, with certain boundary conditions. The grid formation of the plate is done with the meshgrid function.

Afterward, we loop over time t from 0 to tmax with step size dt.

Inside the time loop, we copy the old temperature array T_old to new temperature array T_new for later comparison.

Then, we perform Jacobi Iterations to numerically compute the temperature at each time step. The loop continues until the convergence criterion has been reached, in this case, after 500 iterations.

Finally, we update the T_old with T_new computed and at specific time points, create the contour plots of temperature using the MATLAB contourf function.

Note that the pause(0.1) function at the end of the contour plot inside the for loop is put to slow down the iteration and animation for visualization purposes.

Also, don't forget to define the Jacobi Iterations function.

main.m
% Jacobi's iterations to solve the pressure correction
function jacobi(Tp, Tpo)
    global alpha
    global itermax
    global relresid

    %RHS of the equation
    R=b;

    %Start Jacobi iterations
    for iter = 1:itermax

        %Copy old solution to new array
        Tpo(:,:) = Tp(:,:);

        %Iterate over space (Exclude boundary points)
        for i=2:N
            for j=2:N
                Tp(i,j)=0.25*(Tp(i+1,j)+Tp(i-1,j)+Tp(i,j+1)+Tp(i,j-1))-alpha*Tpo(i,j);
            end
        end

        %Compute the residual and check converge criterion
        resid = norm(R-A*Tp(:));
        if resid < relresid
            break
        end
    end
end
668 chars
30 lines

The above jacobi function computes Jacobi's method with implicit Euler for time, central differences for space. It takes Tp = temperature matrix array, Tpo = temperature matrix array as input and solves for the solution matrix Tp while updating the input Tp array at each iteration.

Overall, this piece of code gives us the contour plot of temperature at several time steps.