set up coefficient matrix and rhs vector using pointer for a 2d unsteady heat conduction equation. implicit euler for time and central differencing for space. use a jacobi solver. in matlab

First, we need to define the problem size and parameters. Let's assume we have a 2D domain of size Nx by Ny, and we want to solve the heat conduction equation with implicit Euler for time and central differencing for space:


main.m
Nx = 50; % number of grid points in x direction
Ny = 50; % number of grid points in y direction
dx = 0.1; % grid spacing in x direction
dy = 0.1; % grid spacing in y direction
k = 1.0; % thermal conductivity
dt = 0.01; % time step
T0 = 0.0; % initial temperature
Tinf = 1.0; % boundary temperature
alpha = k*dt/(dx^2); % spatial diffusion coefficient
351 chars
10 lines

Next, we need to set up the coefficient matrix and the rhs vector using pointers. The coefficient matrix will have (NxNy) rows and columns, and the rhs vector will have (NxNy) elements. To use pointers, we need to define the indices of the matrix elements and rhs vector elements.


main.m
% set up indices of matrix elements and rhs vector elements
ind = @(i,j) (j-1)*Nx + i; % (i,j) -> row index
i0 = @(ind) mod(ind-1,Nx)+1; % index -> i
j0 = @(ind) floor((ind-1)/Nx)+1; % index -> j

% allocate memory for coefficient matrix and rhs vector
A = zeros(Nx*Ny); % coefficient matrix
b = zeros(Nx*Ny,1); % rhs vector

% set up coefficient matrix and rhs vector using central differencing
for j = 1:Ny
    for i = 1:Nx
        % set diagonal element
        A(ind(i,j),ind(i,j)) = 1 + 2*alpha;
        % set left element
        if i > 1
            A(ind(i,j),ind(i-1,j)) = -alpha;
        else % left boundary
            b(ind(i,j)) = b(ind(i,j)) + alpha*Tinf;
        end
        % set right element
        if i < Nx
            A(ind(i,j),ind(i+1,j)) = -alpha;
        else % right boundary
            b(ind(i,j)) = b(ind(i,j)) + alpha*Tinf;
        end
        % set bottom element
        if j > 1
            A(ind(i,j),ind(i,j-1)) = -alpha;
        else % bottom boundary
            b(ind(i,j)) = b(ind(i,j)) + alpha*Tinf;
        end
        % set top element
        if j < Ny
            A(ind(i,j),ind(i,j+1)) = -alpha;
        else % top boundary
            b(ind(i,j)) = b(ind(i,j)) + alpha*Tinf;
        end
    end
end
1247 chars
41 lines

Finally, we can use the Jacobi solver to solve the system of equations. We can define a function for the Jacobi solver, and use it to iterate the solution until convergence.


main.m
% define Jacobi solver function
function u = jacobi(A,b,u0,tol)
% solve the system A*u = b using Jacobi iteration
% A: coefficient matrix (NxN)
% b: rhs vector (Nx1)
% u0: initial guess (Nx1)
% tol: convergence tolerance
N = length(b);
u = u0;
u_new = zeros(N,1);
err = inf;
while err > tol
    % perform one Jacobi iteration
    for i = 1:N
        u_new(i) = (b(i) - A(i,:)*u + A(i,i)*u(i))/A(i,i);
    end
    % calculate error
    err = norm(u_new - u);
    % update solution
    u = u_new;
end
end

% set initial guess for temperature
u0 = T0*ones(Nx*Ny,1);

% iterate the solution until convergence
u = jacobi(A,b,u0,1e-6);
630 chars
29 lines