gistlib - set up coefficient matrix and rhs vector for a 2d unsteady heat conduction equation. implicit euler for time and backward differencing for space in matlab

Assume that we have a 2D unsteady heat conduction problem as follows:


main.m
du/dt = k * (d2u/dx2 + d2u/dy2) + f(x, y, t)
45 chars
2 lines

with initial conditions:


main.m
u(x, y, 0) = u0(x, y)
22 chars
2 lines

and boundary conditions:


main.m
u(x, y, t) = g(x, y, t)    on the boundary of the domain.
58 chars
2 lines

We consider the backward difference method for the spatial derivatives and implicit Euler for time integration with time step size Δt.

Let, Δx and Δy represent the spatial step sizes along x and y directions, respectively. Then, the discretized solution at time nΔt and grid points iΔx and jΔy can be defined as u(i,j,n).

The discrete form of the problem is as follows:


main.m
u(i,j,n+1) = u(i,j,n) + k * Δt / (Δx^2 + Δy^2) * (u(i+1,j,n+1) - 2*u(i,j,n+1) + u(i-1,j,n+1) + u(i,j+1,n+1) - 2*u(i,j,n+1) + u(i,j-1,n+1)) + Δt * f(iΔx,jΔy,nΔt)
161 chars
2 lines

Rearranging the terms, we have:


main.m
-Au(n+1) = Bu(n) + f(n)
24 chars
2 lines

where


main.m
u(n) = [u(1,1,n); u(2,1,n); ...; u(nx,1,n); u(1,2,n); u(2,2,n); ...; u(nx,ny,n)] 

Au(n+1) = [a1,1  a2,1 ...  anx,1  a1,2  a2,2 ...  anx,ny] * [u(1,1,n+1); u(2,1,n+1); ...; u(nx,1,n+1); u(1,2,n+1); u(2,2,n+1); ...; u(nx,ny,n+1)]

B is similar to A with some values set to zero because of boundary conditions.

f(n) is a vector made of f(iΔx,jΔy,nΔt) values at all grid points.    
381 chars
8 lines

The coefficient matrix A and right-hand-side vector B can be generated by performing the following steps:


main.m
% parameters
nx = 10; % number of grid points in x-direction (excluding boundary points)
ny = 10; % number of grid points in y-direction (excluding boundary points)
Lx = 1; % length of domain along x-direction
Ly = 1; % length of domain along y-direction
k = 1; % thermal conductivity
t_end = 10; % simulation end time
dt = 0.01; % time step for simulation
dx = Lx / nx; % grid spacing along x-axis
dy = Ly / ny; % grid spacing along y-axis

% define grid points
x = linspace(0, Lx, nx+2);
y = linspace(0, Ly, ny+2);

% define initial condition
u0 = zeros(ny+2, nx+2);
% TODO: add initial conditions as required

% define boundary conditions
g = zeros(ny+2, nx+2);
% TODO: add boundary conditions as required

% define f(x,y,t) function
f = @(x, y, t) zeros(size(x));

% set up coefficient matrix and rhs vector
I = speye(nx*ny); % identity matrix
A = sparse(nx*ny, nx*ny);
B = sparse(nx*ny, nx*ny);
F = zeros(nx*ny, 1);
for j = 2:ny+1
    for i = 2:nx+1
        p = (j-2)*nx + i-1; % index of current point
        % coefficients for grid point (i,j)
        ax = k*dt/dx^2;
        ay = k*dt/dy^2;
        
        % coefficient for (i,j) in the next time step
        a = 1 + 2*ax + 2*ay;
        
        % coefficient for (i+1,j) in the next time step
        if i < nx+1
            A(p, p+1) = -ax;
        else
            % Dirichlet boundary condition on right boundary
            F(p) = F(p) - ax*g(j,nx+2);
        end
        
        % coefficient for (i-1,j) in the next time step
        if i > 2
            A(p, p-1) = -ax;
        else
            % Dirichlet boundary condition on left boundary
            F(p) = F(p) - ax*g(j,1);
        end
        
        % coefficient for (i,j+1) in the next time step
        if j < ny+1
            A(p, p+nx) = -ay;
        else
            % Dirichlet boundary condition on top boundary
            F(p) = F(p) - ay*g(ny+2,i);
        end
        
        % coefficient for (i,j-1) in the next time step
        if j > 2
            A(p, p-nx) = -ay;
        else
            % Dirichlet boundary condition on bottom boundary
            F(p) = F(p) - ay*g(1,i);
        end
        
        % coefficient for (i,j) in the next time step
        A(p, p) = a;
        
        % coefficient for (i,j) in the current time step
        B(p, p) = 1;
        
        % contribution from source term
        F(p) = F(p) + dt*f(x(i), y(j), 0);
    end
end
2415 chars
84 lines

Finally, we can solve the system of equations for each time step using a linear solver, such as Matlab's backslash operator:


main.m
% initialize solution matrix
U = zeros(ny+2, nx+2);
U(:,1) = g(:,1); % Dirichlet boundary condition on left boundary
U(:,nx+2) = g(:,nx+2); % Dirichlet boundary condition on right boundary
U(1,:) = g(1,:); % Dirichlet boundary condition on bottom boundary
U(ny+2,:) = g(ny+2,:); % Dirichlet boundary condition on top boundary

% loop over time steps
for n = 1:floor(t_end/dt)
    % solve linear system to get u in the next time step
    u_new = reshape(A\(-B*reshape(U(2:ny+1,2:nx+1), [], 1) - F), ny, nx);
    
    % update solution matrix
    U(2:ny+1,2:nx+1) = u_new;
    U(:,1) = g(:,1); % Dirichlet boundary condition on left boundary
    U(:,nx+2) = g(:,nx+2); % Dirichlet boundary condition on right boundary
    U(1,:) = g(1,:); % Dirichlet boundary condition on bottom boundary
    U(ny+2,:) = g(ny+2,:); % Dirichlet boundary condition on top boundary
    
    % TODO: visualize or output solution if required
end
923 chars
22 lines

Note that the above code assumes homogeneous Dirichlet boundary conditions. If you have different boundary conditions, you need to modify the code accordingly. For example, for Neumann boundary conditions, you need to add extra terms to the coefficient matrix and the right-hand-side vector to represent the boundary flux.

set up coefficient matrix and rhs vector for a 2d unsteady heat conduction equation. implicit euler for time and backward differencing for space in matlab

similar matlab code snippets

related categories