measured boundary conditions with pde toolbox
Doug
I see that assempde will take as its input either a boundary condition matrix or a boundary M-file, but both seem to give boundary conditions by specifying a simple formula (like x.^2+y.^2, from the example in the help files). I want to solve Poisson's equation using actual experimental data for the Neumann boundary conditions. It would seem reasonable to provide a matrix similar in form to the edge matrix "e", but containing local values of the Neumann boundary condition. Can it be done?
Thanks in advance from a lowly experimentalist to all you pde gurus!
Bruno Luong
> Has anybody tried using the pde toolbox with boundary conditions taken from experimental measurements (ie, not a formula)?
>
> I see that assempde will take as its input either a boundary condition matrix or a boundary M-file, but both seem to give boundary conditions by specifying a simple formula (like x.^2+y.^2, from the example in the help files). I want to solve Poisson's equation using actual experimental data for the Neumann boundary conditions. It would seem reasonable to provide a matrix similar in form to the edge matrix "e", but containing local values of the Neumann boundary condition. Can it be done?
>
> Thanks in advance from a lowly experimentalist to all you pde gurus!
Unless if you refer to Dirac Neumann bc, the local values Neumann's bc does not make sense. It must be something you can *apply* on another function defined on the boundary. Specifically, mathematician like to refer it as continuous linear application on the fractional Sobolev H^{1/2} space; or even more barbaric: the "H^{-1/2} space". The later does not have defined local value, unfortunately.
Bruno
Doug
Anyway, that's not the sticky part--what I'm unsure about is how to specify boundary conditions from measured data, instead of specifying them as a constant or an analytic function of coordinates.
--Doug
Bruno Luong
> Hmmm, I haven't heard of Dirac-Neumann boundary conditions. In this problem I want to specify the normal derivative of the scalar field at each node on the boundary.
>
> Anyway, that's not the sticky part--what I'm unsure about is how to specify boundary conditions from measured data, instead of specifying them as a constant or an analytic function of coordinates.
>
> --Doug
You can't. Physically, it does not make sense to measure the flux at one point locally. You need to know the flux on the boundary *entirely*. You might assume that your flux is constant or linear by edge, but you measure it at the nodes, but you cannot assume a partial flux. Take a semaphore, you can't estimate accurately the water going out a pond by measuring the water speed at one single point (or many few of them).
Constant Neumann du/dn = g on gamma (the boundary)
<=> specify integral (du/dn.f) dx := integral g.f dx, for all f
Dirac Neuman du/dn = g delta (xi)
<=> specify integral (du/dn.f) dx := g.f(x), for all f
However you won't be able to define the above integral if you only know du/dn at the nodes (you will need to know du/dn in the edge entirely)
Bruno
Doug
Though wbound makes a boundary M-file that is based on a simple formula, in exactly the same format as a boundary condition matrix, it's possible to write your own boundary M-file *any way you like*. As long as it takes (p,e,u,time) as inputs and returns [q,g,h,r] as outputs, it works fine with assempde. You just pass the filename of the boundary M-file as the first argument to assempde. So I've written a boundary M-file that takes the necessary inputs (p,e,u,time) but ignores them, instead constructing [q,g,h,r] from my measurements. And it works. Reading the help file for pdebound over and over was the key!
Bruno, I'm still not following your argument that Neumann boundary conditions can't work. I understand that they leave room for an arbitrary constant of integration, but that's OK for me, because in my application, the solution to Poisson's equation is a stream function. My real interest is its gradient, for which the constant of integration is irrelevant. You seem to be concerned about the fact that the nodes are discrete points, but so are my measurements, and so is every numerical grid; again I don't understand. In case we're still talking about two different things, here's what I mean by Neumann boundary conditions:
http://en.wikipedia.org/wiki/Neumann_boundary_condition
http://mathworld.wolfram.com/BoundaryConditions.html
Bruno Luong
You never get the right convergence property with *point-wise* Neumann condition, i.e., your PDE is ill posed; from that you might get all nasty stuff, and the numerical solution might be very far from the true solution, whatever it is.
But nevermind, nobody never get the sense why we bother to define all these Sobolev space anyway, it is not just for fun.
Bruno
Bruno Luong
I might not fully understand where your Neumann data are from (I understand from measurement), but you have to be careful about how to plug them into the PDE, even if you think you can define it in the boundary.m.
Bruno
Carl S.
>
> I might not fully understand where your Neumann data are from (I understand from measurement), but you have to be careful about how to plug them into the PDE, even if you think you can define it in the boundary.m.
>
> Bruno
u=NeumannBoundCond(u);
[u_x,u_y]=gradient(u);
s=sqrt(u_x.^2 + u_y.^2);
smallNumber=1e-10;
Nx=u_x./(s+smallNumber); % add a small positive number to avoid division by zero
Ny=u_y./(s+smallNumber);
curvature=div(Nx,Ny);
diracPhi=Dirac(u,epsilon);
areaTerm=diracPhi.*g; % balloon/pressure force
edgeTerm=diracPhi.*(vx.*Nx+vy.*Ny) + diracPhi.*g.*curvature;
% u=u + timestep*(lambda*EdgeTerm + alfa*areaTerm);
SobolevGradEdgeTerm = (1./(Img-gradient(u).^2)).*edgeTerm;
u=u + timestep*(lambda*SobolevGradEdgeTerm + alfa*areaTerm);
When I multiply the edgeTerm and lambda then I can see the evoluation of the active contour. But I want to use Sobolev gradient. (So, I have added the last two code lines.) In this case, when I multiply the lambda with SobolevGradEdgeTerm then the active contour disappears in the figure. What is the mistake ?
Carl S.
SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
but I have still the same problem :(
Carl S.
1./(eye(256)-gradient(u).^2)) becomes inf. Therefore, the active contour disappears
How to solve this problem ?
Carl S.
Carl S.
temp2 = 1./(eye(size(u))-gradient(u).^2);
temp2(isequal(temp2, inf)) = 1e-10;
but temp2 have still inf values, why ? please help me
Carl S.
[row,col] = find(temp2==Inf);
temp2(row,col)=1e-10;
but, there is oversegmentation problem :(((((((
ImageAnalyst, where are you ? Any idea, any sugesstion ?
Dmitry
>
> Though wbound makes a boundary M-file that is based on a simple formula, in exactly the same format as a boundary condition matrix, it's possible to write your own boundary M-file *any way you like*. As long as it takes (p,e,u,time) as inputs and returns [q,g,h,r] as outputs, it works fine with assempde. You just pass the filename of the boundary M-file as the first argument to assempde. So I've written a boundary M-file that takes the necessary inputs (p,e,u,time) but ignores them, instead constructing [q,g,h,r] from my measurements. And it works. Reading the help file for pdebound over and over was the key!
Could you provide an example of this kind of usage please?
As I understood, you just loaded the variable in the boundary M-file like this:
> load('var.mat','var');
and used it for construction of [q,g,h,r]. Am I right?
I've done the same (for hyperbolic solver), but I can't check if it works because for some (another) reason the solver gives me the error:
>Error using pdeexpd (line 57)
>Number of variables mismatch between u and bl.
>Error in assemb (line 101)
> [q,g,h,r]=pdeexpd(p,e,u,time,bl);
>Error in pdeODEInfo/getMats (line 164)
> [Q,G,H,R]=assemb(self.b,self.p,self.e,u,time);
>Error in pdeODEInfo/checkFuncDepen (line 60)
> [MM0,K0,M0,F0,Q0,G0,H0,R0] = self.getMats(u0, t0);
>Error in pdeHyperbolicInfo (line 15)
> obj=obj.checkFuncDepen(u0, tlist);
>Error in hyperbolic (line 41)
> pdehyp=pdeHyperbolicInfo(u0,ut0,tlist,b,p,e,t,c,a,f,d);
>Error in main (line 49)
>u = hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);
In my case I have only dependence on time in boundary M-file, but the solver assumes that output of it depends on the solution u too
(it goes wrong way at pdeODEInfo/checkFuncDepen stage).
Don't know what to do...
Bill Greene
with how you are passing the handle to it into hyperbolic.
Did you follow the example here if you have a single pde:
http://www.mathworks.com/help/pde/ug/boundary-conditions-for-scalar-pde.html
or here if you have a system of pde:
http://www.mathworks.com/help/pde/ug/boundary-conditions-for-pde-systems.html
Is the value of "b" that you are passing into hyperbolic something like this?
b = @myPdeBoundFunction; % myPdeBoundFunction is whatever you named your function
Bill
Dmitry
Dmitry
I have a time-dependent boundary conditions for hyperbolic solver, but when I displayed the variable time, which is passed in boundary.m, I saw it was nearly equal to 2*10^-15 during all calls. So hyperbolic() passes time=0 to BC-function all the time. Which it can be related to?
Dmitry
I realized that I need to add
> if isnan(time)
> gmatrix = ones(1,ne)*NaN;
> end
for make pdebound-file time dependent. Now everything is ok.