Specrtal Method for Heat Equation
1 次查看(过去 30 天)
显示 更早的评论
I have one dimensional homogeneous heat eqaution.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248968/image.png)
I want to solve it numerically using by supposing
where ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248970/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248969/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248970/image.png)
Taking time derivative using farward difference method
, after combinging I get
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248971/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248972/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248973/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248974/image.png)
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 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Subplots 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!