Solving the heat equation using finite difference with tridiagonal matrix
12 次查看(过去 30 天)
显示 更早的评论
I have the following code to solve the heat equation via a tridiagonal matrix but have no clue how to continue.
clear
nx=10;
a=0; b=1;
Du=1; dt = 0.1;
alpha=1; beta=1; %BCs
p=0; q=0;
x=linspace(a,b,nx+1); dx=(b-a)/nx; %forms the grid and computes the step-size
u(1)=alpha; u(nx+1)=beta; %sets the boundary conditions
A=tridiag(1/(dx*dx)-p/(2*dx),q-2/(dx*dx),1/(dx*dx)+p/(2*dx),nx-1);
% A is the tri-diagonal matrix containing the LHS of the finite diff equations
r=functionr(x)';
% r gives the RHS function.
% The Dirichlet bc must be incorporated in the 1st and
% last equation of the system as follows:
r(2)=r(2)-(1/(dx*dx)-p/(2*dx))*alpha;
r(nx)=r(nx)-(1/(dx*dx)+p/(2*dx))*beta;
% u(1) and u(n+1) given by BC,
% values in between computed from solving the linear system as follows:
% Euler's Method
h = (Du*dt)/dx^2;
for j=1:nx-1 % loop for Euler's method
u(j+1) = u(j) + h*A*u(j);
end
plot(x,u)
the tridiag function is as follows
function T = tridiag(a, b, c, n)
% tridiag Tridiagonal matrix.
% T = tridiag(a, b, c, n) returns an n by n matrix that has
% a, b, c as the subdiagonal, main diagonal, and superdiagonal
% entries in the matrix.
T = b*diag(ones(n,1))+c*diag(ones(n-1,1),1)+a*diag(ones(n-1,1),-1);
the functionr is as follows
function r = functionr(x)
r = x*0;
all help much appreciated
回答(1 个)
Rohit Kulkarni
2024-3-20
Hi Joshua,
In my understaing it looks like you're trying to solve a differential equation using a finite difference method, but there seems to be some issues in your code:
Issues and Improvements:
1. Initialization of `u`:You need to initialize `u` for all grid points, not just the boundaries. Initially, you can set it to zeros or any initial condition if given.
2. Matrix-vector multiplication: In the Euler method loop, the way you're trying to update `u(j+1)` is incorrect. You should perform a matrix-vector multiplication for the whole vector `u` (excluding the boundary values) at once, not element by element.
3. Boundary conditions handling: The boundary conditions are correctly set at the beginning, but when updating `u`, you should exclude the boundaries since they are constant (for Dirichlet conditions).
4. Function functionr:The function `functionr` returns a zero vector, which means it doesn't contribute to the solution. If there's supposed to be a source term or any non-zero right-hand side, it should be defined here.
Improved Code:
clear
nx = 10;
a = 0; b = 1;
Du = 1; dt = 0.1;
alpha = 1; beta = 1; % BCs
p = 0; q = 0;
x = linspace(a, b, nx+1); dx = (b-a)/nx; % forms the grid and computes the step-size
u = zeros(nx+1,1); % Initialize u for all grid points
u(1) = alpha; u(nx+1) = beta; % sets the boundary conditions
% Define A matrix
A = tridiag(1/(dx*dx)-p/(2*dx), q-2/(dx*dx), 1/(dx*dx)+p/(2*dx), nx-1);
% Define r vector
r = functionr(x(2:end-1))'; % Exclude the boundary points in x
% Adjust r for boundary conditions
r(1) = r(1) - (1/(dx*dx)-p/(2*dx))*alpha;
r(end) = r(end) - (1/(dx*dx)+p/(2*dx))*beta;
% Time-stepping with Euler's method
h = (Du*dt)/dx^2;
for j = 1:nx-1
u(2:nx) = u(2:nx) + h*(A*u(2:nx) + r); % Update interior points only
end
plot(x, u)
xlabel('x')
ylabel('u')
title('Solution of the Heat Equation')
Hope this resolves your query!
Thanks
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!