Merge branch 'master' of https://git.art-ist.cc/daniel/NSSC
This commit is contained in:
commit
a4bed684f7
|
|
@ -0,0 +1,59 @@
|
|||
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
|
||||
Binary file not shown.
|
|
@ -0,0 +1,63 @@
|
|||
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
|
||||
|
|
@ -0,0 +1,62 @@
|
|||
function [x,t,u] = Mixed_EI(L,N,T,K,c1,c2,f,u0,D)
|
||||
% ---- Risoluzione dell'equazione del calore ----
|
||||
% u_t - u_xx = f nell'intervallo [-L,L]
|
||||
% con condizioni al bordo di Dirichlet
|
||||
% e condizioni iniziali.
|
||||
% -----------------------------------------------
|
||||
% Sintassi:
|
||||
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
|
||||
%
|
||||
% Input:
|
||||
% L semiampiezza intervallo spaziale (-L,L)
|
||||
% N numero di sottointervalli in (-L,L)
|
||||
% T estremo finale intervallo temporale (0,T)
|
||||
% K numero di sottointervalli in (0,T)
|
||||
% c1 funzione che descrive la condizione di Dirichlet in x=-L
|
||||
% c2 funzione che descrive la condizione di Dirichlet in x=L
|
||||
% f funzione che descrive il termine noto dell'equazione
|
||||
% u0 funzione che descrive la condizione iniziale in t=0
|
||||
%
|
||||
% Output:
|
||||
% x vettore dei nodi spaziali
|
||||
% t vettore dei nodi temporali
|
||||
% u soluzione numerica
|
||||
% del problema
|
||||
|
||||
% Calcolo passo di discretizzazione in spazio e tempo
|
||||
h=L/N;
|
||||
tau=T/K;
|
||||
% Inizializzazione del vettore t
|
||||
t=linspace(0,T,K+1);
|
||||
% Inizializzazione del vettore x
|
||||
x=linspace(0,L,N+1);
|
||||
|
||||
% Inizializzazione della matrice soluzione u
|
||||
u=zeros(N+1,K+1);
|
||||
% Condizione iniziale
|
||||
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);
|
||||
|
||||
|
||||
for k=1:K
|
||||
% Assemblaggio termine noto
|
||||
F=f(x(2:end-1),t(k+1));
|
||||
% Correzione del termine noto con le condizioni al bordo
|
||||
F(1)=F(1) + c1(t(k+1))/(h^2);
|
||||
F(end)=c2*h*2;
|
||||
% Risoluzione del problema
|
||||
u(2:end-1,k+1) = ((I+tau*A)\u(2:end-1,k))' + tau*F;
|
||||
u(end,k+1)=u(end-2,k+1) + c2*2*h;
|
||||
end
|
||||
|
|
@ -0,0 +1,93 @@
|
|||
%% eulero forward DIrichlet
|
||||
L=pi;
|
||||
T=2;
|
||||
f=@(x,t) 0*x.*t;
|
||||
c1=@(t) 1+0*t;
|
||||
c2=@(t) 0*t;
|
||||
u0=@(x) 0*x;
|
||||
D=0.5;
|
||||
%uex=@(x,t) cos(x).*exp(t);
|
||||
|
||||
N=10;
|
||||
K=200;
|
||||
[x,t,u]=Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D);
|
||||
|
||||
figure(1)
|
||||
for ii=1:K+1
|
||||
plot(x,u(:,ii)');
|
||||
xlim([0 L])
|
||||
pause(0.05);
|
||||
end
|
||||
|
||||
space=linspace(0,L,101);
|
||||
time=linspace(0,T,201);
|
||||
[xx,yy]=meshgrid(time,space);
|
||||
%exsol=uex(yy,xx);
|
||||
figure(2)
|
||||
mesh(t,x,u)
|
||||
%figure(2)
|
||||
%mesh(xx,yy,exsol)
|
||||
|
||||
%% Eulero forward Mixed BC
|
||||
|
||||
L=2*pi;
|
||||
T=5;
|
||||
f=@(x,t) 0*x.*t;
|
||||
c1=@(t) 1+0*t;
|
||||
c2=0;
|
||||
u0=@(x) 0*x;
|
||||
D=1.1;
|
||||
%uex=@(x,t) cos(x).*exp(t);
|
||||
|
||||
N=25;
|
||||
K=200;
|
||||
[x,t,u]=Mixed_EA(L,N,T,K,c1,c2,f,u0,D);
|
||||
|
||||
figure(1)
|
||||
for ii=1:K+1
|
||||
plot(x,u(:,ii)');
|
||||
xlim([0 L])
|
||||
ylim([0 1.5])
|
||||
pause(0.02);
|
||||
end
|
||||
|
||||
space=linspace(0,L,101);
|
||||
time=linspace(0,T,201);
|
||||
[xx,yy]=meshgrid(time,space);
|
||||
%exsol=uex(yy,xx);
|
||||
figure(2)
|
||||
mesh(t,x,u)
|
||||
%figure(2)
|
||||
%mesh(xx,yy,exsol)
|
||||
|
||||
%% Eulero Backward Mixed BC
|
||||
|
||||
L=2*pi;
|
||||
T=5;
|
||||
f=@(x,t) 0*x.*t;
|
||||
c1=@(t) 1+0*t;
|
||||
c2=0;
|
||||
u0=@(x) 0*x;
|
||||
D=1;
|
||||
%uex=@(x,t) cos(x).*exp(t);
|
||||
|
||||
N=25;
|
||||
K=200;
|
||||
[x,t,u]=Mixed_EI(L,N,T,K,c1,c2,f,u0,D);
|
||||
|
||||
figure(1)
|
||||
for ii=1:K+1
|
||||
plot(x,u(:,ii)');
|
||||
xlim([0 L])
|
||||
ylim([0 1.5])
|
||||
pause(0.02);
|
||||
end
|
||||
|
||||
space=linspace(0,L,101);
|
||||
time=linspace(0,T,201);
|
||||
[xx,yy]=meshgrid(time,space);
|
||||
%exsol=uex(yy,xx);
|
||||
figure(2)
|
||||
mesh(t,x,u)
|
||||
%figure(2)
|
||||
%mesh(xx,yy,exsol)
|
||||
Loading…
Reference in New Issue