NSSC/Exercise_03/Dirichlet_EA.m

60 lines
1.5 KiB
Matlab
Executable File

function [x,t,u] = Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D)
% ---- Solution of heat equation ----
% u_t - u_xx = f in the interval [-L,L] (doesn't matter if we change it)
% Dirichlet BC
%
% -----------------------------------------------
% Sintassi:
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
%
% Input:
% L Half of the width (-L,L)
% N number of intervals in (-L,L)
% T max time (0,T)
% K number of intervals in (0,T)
% c1 Dirichlet BC in x=-L
% c2 Dirichlet BC in x=L
% f force function
% u0 initial condition in t=0
%
% Output:
% x vector of the spatial nodes
% t vector of the time nodes
% u numeric solution of the problem
% discretisation step in time and space
h=L/N; %space
tau=T/K; %time
% Initialization of t
t=linspace(0,T,K+1)';
% initialization of x
x=linspace(0,L,N+1);
% Initialization of the matrix solution u
u=zeros(N+1,K+1);
% Intial conditions
u(:,1)=u0(x);
% BC
u(1,:)=c1(t);
u(end,:)=c2(t);
% Creation of the matrix A
% We obtained this matrix using the Central Discretization method
e=ones(N-1,1);
A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
I=speye(N-1,N-1);
% Compute the solution
% What we are doing here is to compute the solution in the interval for
% each time k
for k=1:K
% Assembly of the force
F=f(x(2:end-1),t(k));
% Correction of the force using the BC
F(1)=F(1) + D*c1(t(k))/(h^2);
F(end)=F(end) + D*c2(t(k))/(h^2);
% Solution using Eulero forward
u(2:end-1,k+1) = ((I - tau*D*A)*u(2:end-1,k))' + tau*F;
end