pdepe problem (numerical vs analytical solution)
3 次查看(过去 30 天)
显示 更早的评论
So this is my pde:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/475873/image.png)
The IC for 0<x<1:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/476398/image.png)
The b.c for t>0.
u is bounded.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/475883/image.png)
Using separation of variables by letting v = u-1, we then have the analytical solution:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/476403/image.png)
The code in pdepe:
function [c,f,s] = pdex1pde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = 0; %ignored because m=2
ql = 0; %ignored because m = 2
pr = ur-1;
qr = 0;
end
function u0 = u0(x)
u0 = 0; %initial condition
end
m = 2;
sol = pdepe(m,@pdex1pde,@u0,@pdex1bc,x,t);
u = sol(:,:,1);
x = linspace(0,1,50); %x = linspace(0,L,mesh points);
t = linspace(0,2,50);
Firstly, I do not see exponential decay or sinusoidal waves that should be expected from the analytical solution.
surf(x,t,u)
title('Numerical solution computed with 50 mesh points')
xlabel('Distance x')
ylabel('Time t')
Also, my analytical vs numerical solution does not match up.
Could someone tell me how to fix this?
plot(x,u(25,:),'o',x,1-(2*x.^(-1)*(pi)^(-1))'*exp(-pi^2*t(25))'*sin(pi*x))
title('Solution at t = 1')
legend('Numerical, 50 mesh points','Analytical','Location','South')
xlabel('Distance x')
ylabel('u(x,1)')
0 个评论
采纳的回答
Pat Gipper
2021-1-3
I don't think you were summing the terms of your series expansion correctly. The attached code carries out the expansion summation to 900 terms which seemed adequate. I also changed the source coefficient expression to what I thought is the correct equation. But regardless whether the source is your original, "s=0", or my answer, "s=(2/x) * dudx", there is still a mismatch between the series expansion and numerical solution. Oddly, I found that if you used "s=(1/x) * dudx" for the source there is very little error. Additionally, the results were no different whether a symmetry constant of m=1 or m=2 was selected.
In summary, I think the numerical solution is giving results that are making sense. We simply need to decide what source coefficient expression should be selected.
0 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!