pdepe problem (numerical vs analytical solution)

조회 수: 2 (최근 30일)
Ricardo Machado
Ricardo Machado 2021년 1월 1일
답변: Pat Gipper 2021년 1월 3일
So this is my pde:
The IC for 0<x<1:
The b.c for t>0.
u is bounded.
Using separation of variables by letting v = u-1, we then have the analytical solution:
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)')

채택된 답변

Pat Gipper
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개)

카테고리

Help CenterFile Exchange에서 Eigenvalue Problems에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by