Specrtal Method for Heat Equation
    6 views (last 30 days)
  
       Show older comments
    
I have one dimensional homogeneous heat eqaution.

 I want to solve it numerically using by supposing  where
where 
 where
where 
Taking time derivative using farward difference method  , after combinging I get
 , after combinging I get
 , after combinging I get
 , after combinging I get , where
, where  and A is sparse matrix with
 and A is sparse matrix with  on its main diagonal.
 on its main diagonal.I tried to write it in MATLAB, but got wrong answers, please look into my code
clc;clear all;close all;
%  problem    u_t =  beta u_xx%     beta 1
%  exact solution
%%
u  = @(x,t) (4*pi/3)*(((sin(pi*x))*exp(-(pi^2)*t)));
t0  = 0;
tn  = 0.8;
x0  = 0;
xn  = 1;
dx = 0.0625;
dt = 0.004;
x = (x0:dx:xn)';
t = (t0:dt:tn)';
beta =  1;
s    =  beta*dt/dx^2;
nx = (xn-x0)/dx;
nt = (tn-t0)/dt;
%% we start indexing from zero
nx_int  = nx;  %  number of interior points in spatial dim
nt_int = nt-1;  %  number of interior points in temporal dim
%% construction of diagonal matrix
indx1 = ones(1,nx-2);
indx2 = ones(1,nx-1);
v =((1:nx)*pi).^2;
A = full(diag(v,0));
%% boundary conditions
%  U(0,k)  = 0
%  U(20,k) = 0
% initial condition
% int_cond = @(x) sin(pi/4*x).*(1+2*cos(pi/4*x));
%%
fun = @(x) x.*(1-x).*sin((1:nx+1)*pi*x);
q = integral(fun,0,1,'ArrayValued',true);
U=zeros(nx+1,nt+1);
U_exact=zeros(nx+1,nt+1);
U(:,1) = q';
U(1,:)=0;
U(nx+1,:)=0;
% A U(k+1) = U(k) + b
%% numerical solution construction
for k=1:nt_int+1
    U(2:nx_int+1,k+1) = U(2:nx_int+1,k)+ dt*A\U(2:nx_int+1,k+1);
end
%% exact solution construction
for k=1:nt+1
    U_exact(1:nx+1,k) =u(x,t(k));
end
Numerical_sol = U;
Analytical_sol = U_exact;
%Plots
figure(1)
subplot(2,2,1)
mesh(t,x,U)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Numerical sol')
subplot(2,2,2)
mesh(t,x,U_exact)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Analytical sol')
subplot(2,2,3)
mesh(t,x,U_exact-U)
xlabel('Temporal dim')
ylabel('Spatial dim')
zlabel('Error in  sol')
and tell me my mistake.
Thanks
0 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!