Specrtal Method for Heat Equation

1 visualización (últimos 30 días)
Muhammad Usman
Muhammad Usman el 19 de Nov. de 2019
I have one dimensional homogeneous heat eqaution.
I want to solve it numerically using by supposing where
Taking time derivative using farward difference method , after combinging I get
, where and A is sparse matrix with 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

Respuestas (0)

Categorías

Más información sobre Surface and Mesh Plots en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by