Solving an elliptic PDE with a point source
4 views (last 30 days)
Show older comments
How do I solve the following PDE:
where c is a constant and delta is a point source? The domain is a rectangle and the boundary conditions are that u is zero on the boundary and the second derivative of u normal to the boundary is also zero. How do I cast this equation in terms of the generic MATLAB PDE,
which is solved by u = assempde(b,p,e,t,c,a,f)?
3 Comments
Bruno Pop-Stefanov
on 21 Jan 2014
I think you can follow a procedure similar to that example for the heat equation:
Give me a little more time to try it myself.
Answers (2)
Youssef Khmou
on 21 Jan 2014
Edited: Youssef Khmou
on 21 Jan 2014
Alex, You can find many tutorials on how to use the PDETOOL, however i recommend that you solve the equation by program , later you can verify your example, i tried to write the finite differences method for your equation, try to manipulate it a little a bit :
% PDE
clear;
N=200;
U=zeros(N);
Delta=1;
U(N/2,N/2)=Delta;
c=2;
T=400;
for t=1:T
for x=2:N-1
for y=2:N-1
if U(x,y)==Delta
continue;
else
E=(U(x+1,y)+U(x-1,y)+U(x,y+1)+U(x,y-1));
U(x,y)=((Delta/c)-E-((1/c)*U(x+1,y)))/(-(4+(1/c)));
end
end
end
end
figure, contour(U), shading interp
2 Comments
Youssef Khmou
on 22 Jan 2014
Alex , t is considered as the number of iterations necessary for convergence such |U(x,y)t+1|-||U(x,y)t||<tolerance, anyway we wait for an answer using pdetool .
Cordially .
Bill Greene
on 22 Jan 2014
Hi,
As you've probably observed, PDE Toolbox doesn't support the du/dx term explicitly. However there is a trick you can try that often works.
First, I suggest you start with this example:
This example shows how to approximate the delta function and use adaptive meshing to refine the mesh around the load.
The trick to dealing with the du/dx term is to move it to the RHS and define the f-coefficient to be a function of ux (aka du/dx).
I made the following changes to the example I mention above:
I added these options to the call to adaptmesh:
'Nonlin', 'on', 'jac', 'full'
Then I made these changes to the circlef function:
Added:
[ux,uy] = pdegrad(p, t, u);
f = -ux;
and changed
if ~isnan(tn)
f(tn) = f(tn) + 1/ar(tn);
end
Bill
0 Comments
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!