Convergence of Heat Equation?
조회 수: 1 (최근 30일)
이전 댓글 표시
%% Ut=Uxx,
%%solution: U(x,t)=sin(2*pi*x)*exp(-(2*pi).^2*t),
%%I.C : U(x,0)= sin(2*pi*x)
%%B.C: U(0,t)=U(1,t)=0,
L = 1 ; %Length of Domain
t_final = 0.1; % time limit
N = 32; % no of elements
dx = 1/N; % step size
dt = 0.001
% x = 0:dx:L; %spatial discretization
x = linspace(0,L,N+1);
t = 0:dt:t_final; %temporal discretization
a= length(x);
b = length(t);
U = zeros(1,b);
W = zeros(a,b);
X = sin(2*pi*x);
W(:,1) =X;
W(:, b)=0;
% t=0;
T = zeros(b,1);
T(1) = 0;
for i=1:length(t)
T= T+dt;
U = sin(2*pi*x)*(exp((-(2*pi).^2)*T(i,:)));
W(:,i) = U;
X = U;
U(1) =0;
U(end)=0;
T(i) = T(1);
figure(1)
plot(x,U)
hold on
end
%% Backward Euler for FEM solution
xx = linspace(0,L,N+1); %Fem spatial discretization
% xx = 0:dx:L;
nx = length(xx); %total nodes
M = spdiags([(dx/6)*ones(nx,1), (2*dx/3)*ones(nx,1), (dx/6)*ones(nx,1)],[-1, 0,1],nx,nx) ; %% Mass matrix
A = spdiags([(-1/dx)*ones(nx,1), (2/dx)*ones(nx,1), (-1/dx)*ones(nx,1)],[-1, 0, 1],nx,nx) ;%% Stiffness matrix
delta_t = 0.001;
m = t_final/delta_t;
delta_t = t_final/m;
tspan = 0:delta_t:t_final;
% store solution for each time in matrix u
u = zeros(nx,m+1);
tvec = zeros(1,m+1);
U_tk = sin(2*pi*xx);
u(:,1) = U_tk;
tt = 0 ;
tvec(1) = tt;
tt = 0;
b = zeros(1,nx); % Neumann boundary condition
b(nx) = 0;
for nt = 1:m
tt = tt+delta_t;
c = U_tk + b;
w = (M+delta_t*A)\M;
FEM = c*w;
FEM(1) = 0;
FEM(end)= 0;
u(:, nt) = FEM;
% for next time step:
U_tk = FEM;
tvec(nt) = tt;
figure(2)
plot(xx,FEM)
end
plot(x,U-FEM); )
ERROR = sqrt(sum( dx*(U-FEM).^2 ))
댓글 수: 0
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Heat and Mass Transfer에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!