64 lines
1.7 KiB
Mathematica
64 lines
1.7 KiB
Mathematica
|
function [x,t,u] = Mixed_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 Dirichlet
|
||
|
|
u(1,:)=c1(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);
|
||
|
|
%modify the matrix such that we can use the noimann condition to calculate
|
||
|
|
%the last node.
|
||
|
|
%In this case we want to use a second order decentralized approximation
|
||
|
|
A(end,end-1)=-2/(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) + 2*h*c2;
|
||
|
|
% Solution using Eulero forward
|
||
|
|
u(2:end-1,k+1) = ((I - D*tau*A)*u(2:end-1,k))' + tau*F;
|
||
|
|
u(end,k+1)=u(end-2,k+1) + c2*2*h*D;
|
||
|
|
end
|