2014-09-24 52 views
0

我试图解决1D中的瞬态热方程,并比较分析解和数值解。这些解决方案具有相同的趋势但非常不同,即使其明显不应该是相对误差也是零。我不知道我是否已经正确解决了PDE。 (pde是du/dt = d^2u/dx^2)并且bcs是u(0,t)= 1,u(100,t)= 0,u(x,0)= 0。有人可以看看我的代码吗?Matlab PDE求解器问题

function he 
m = 0; 
x = linspace(0,100,500); 
t = linspace(0,1000,500); 
sol = pdepe(m,@hepde,@heic,@hebc,x,t); 
u = sol(:,:,1); 
y = erfc(x./(2*(t.^0.5))); 
r=(y-u(70,:))/y; 
figure; 
plot(x,u(50,:),'.',x,u(150,:),'.',x,u(250,:),'.',x,u(end,:),'.',x,y,'.'); 
title('Numerical Solutions at different times.'); 
legend('t=100','t=300','t=500','t=700','y ana',0); 
xlabel('Distance x'); 
ylabel('u(x,t)'); 
figure; 
plot(x,r); 
title('error in numerical and analytical solution'); 
legend('error',0); 
xlabel('Distance x'); 
ylabel('error'); 

% -------------------------------------------------------------------------- 
function [c,f,s] = hepde(x,t,u,DuDx) 
c = 1; 
f = DuDx; 
s = 0; 
% -------------------------------------------------------------------------- 
function u0 = heic(x) 
u0 = 0; 
% -------------------------------------------------------------------------- 
function [pl,ql,pr,qr] = hebc(xl,ul,xr,ur,t) 
pl = ul-1; 
ql = 0; 
pr = ur; 
qr = 0; 

回答

0

实际上您的解决方案都是正确设置的,您只是错误地使用了结果。我会一一浏览错误。

sol = pdepe(m,@hepde,@heic,@hebc,x,t); 
u = sol(:,:,1); 

你找到正确的答案,但线路u=sol(:,:,1);是无用的size(sol)=[2 2],所以你可能也只是做u=pdepe(...);

现在,当你计算出确切的解决方案时,你做得非常奇怪。你想在每个x/t的组合中找到它,但你只在其中的一部分组合。您需要使用meshgrid以获得xt的所有组合,然后计算每个组合的确切解。

[X,T]=meshgrid(x,t); 
y = erfc(X./(2*(T.^0.5))); 

然后您需要以不同的方式计算错误。

r=(y-u)./y; 
figure(3); 

而且情节不同。

plot(x,u(50,:),'b',x,u(150,:),'g',x,u(250,:),'r',x,u(end,:),'c',x,y(50,:),'b--',x,y(150,:),'g--',x,y(250,:),'r--',x,y(end,:),'c--'); 

其实你确切的解决方案不会在x=100满足边界条件......

+0

太感谢你了大卫,我得到了解决,现在!!!! – 2014-09-25 03:01:13