Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
measured boundary conditions with pde toolbox

Subject: measured boundary conditions with pde toolbox

From: Doug

Date: 5 Aug, 2009 19:36:04

Message: 1 of 18

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!

Subject: measured boundary conditions with pde toolbox

From: Bruno Luong

Date: 5 Aug, 2009 22:21:01

Message: 2 of 18

"Doug " <dhk@umd.edu> wrote in message <h5cmv4$lbc$1@fred.mathworks.com>...
> 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

Subject: measured boundary conditions with pde toolbox

From: Doug

Date: 6 Aug, 2009 01:16:08

Message: 3 of 18

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

Subject: measured boundary conditions with pde toolbox

From: Bruno Luong

Date: 6 Aug, 2009 05:27:03

Message: 4 of 18

"Doug " <dhk@umd.edu> wrote in message <h5daso$lqr$1@fred.mathworks.com>...
> 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

Subject: measured boundary conditions with pde toolbox

From: Doug

Date: 10 Aug, 2009 19:36:19

Message: 5 of 18

OK, I've answered my own question, and the answer is a boundary M-file.

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

Subject: measured boundary conditions with pde toolbox

From: Bruno Luong

Date: 10 Aug, 2009 19:54:04

Message: 6 of 18

"Doug " <dhk@umd.edu> wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> OK, I've answered my own question, and the answer is a boundary M-file.
>
> 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

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

Subject: measured boundary conditions with pde toolbox

From: Bruno Luong

Date: 11 Aug, 2009 09:52:05

Message: 7 of 18

Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.

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

Subject: measured boundary conditions with pde toolbox

From: Carl S.

Date: 17 Nov, 2012 06:14:09

Message: 8 of 18

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
>
> 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

Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;

    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 ?

Subject: measured boundary conditions with pde toolbox

From: Carl S.

Date: 17 Nov, 2012 14:02:16

Message: 9 of 18

"Carl S." wrote in message <k879vh$k4n$1@newscl01ah.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> > Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
> >
> > 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
>
> Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;
>
> 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 ?

I tried this
SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
but I have still the same problem :(

Subject: measured boundary conditions

From: Carl S.

Date: 18 Nov, 2012 04:37:19

Message: 10 of 18

"Carl S." wrote in message <k885d7$h85$1@newscl01ah.mathworks.com>...
> "Carl S." wrote in message <k879vh$k4n$1@newscl01ah.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> > > Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
> > >
> > > 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
> >
> > Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;
> >
> > 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 ?
>
> I tried this
> SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
> but I have still the same problem :(

(eye(256)-gradient(u).^2) becomes zero. Therefore,
1./(eye(256)-gradient(u).^2)) becomes inf. Therefore, the active contour disappears

How to solve this problem ?

Subject: measured boundary conditions

From: Carl S.

Date: 18 Nov, 2012 07:29:12

Message: 11 of 18

"Carl S." wrote in message <k89olv$rr$1@newscl01ah.mathworks.com>...
> "Carl S." wrote in message <k885d7$h85$1@newscl01ah.mathworks.com>...
> > "Carl S." wrote in message <k879vh$k4n$1@newscl01ah.mathworks.com>...
> > > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> > > > Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
> > > >
> > > > 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
> > >
> > > Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;
> > >
> > > 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 ?
> >
> > I tried this
> > SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
> > but I have still the same problem :(
>
> (eye(256)-gradient(u).^2) becomes zero. Therefore,
> 1./(eye(256)-gradient(u).^2)) becomes inf. Therefore, the active contour disappears
>
> How to solve this problem ?

If I assign a small value that is 1e-10.*randi(1,1) instead of zeros then the above codes causes oversegmentation :(

Subject: measured boundary conditions

From: Carl S.

Date: 18 Nov, 2012 08:26:13

Message: 12 of 18

"Carl S." wrote in message <k8a2o8$2sj$1@newscl01ah.mathworks.com>...
> "Carl S." wrote in message <k89olv$rr$1@newscl01ah.mathworks.com>...
> > "Carl S." wrote in message <k885d7$h85$1@newscl01ah.mathworks.com>...
> > > "Carl S." wrote in message <k879vh$k4n$1@newscl01ah.mathworks.com>...
> > > > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> > > > > Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
> > > > >
> > > > > 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
> > > >
> > > > Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;
> > > >
> > > > 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 ?
> > >
> > > I tried this
> > > SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
> > > but I have still the same problem :(
> >
> > (eye(256)-gradient(u).^2) becomes zero. Therefore,
> > 1./(eye(256)-gradient(u).^2)) becomes inf. Therefore, the active contour disappears
> >
> > How to solve this problem ?
>
> If I assign a small value that is 1e-10.*randi(1,1) instead of zeros then the above codes causes oversegmentation :(

I have tried this
      temp2 = 1./(eye(size(u))-gradient(u).^2);
      temp2(isequal(temp2, inf)) = 1e-10;
but temp2 have still inf values, why ? please help me

Subject: measured boundary conditions

From: Carl S.

Date: 18 Nov, 2012 11:36:08

Message: 13 of 18

"Carl S." wrote in message <k8a635$dg7$1@newscl01ah.mathworks.com>...
> "Carl S." wrote in message <k8a2o8$2sj$1@newscl01ah.mathworks.com>...
> > "Carl S." wrote in message <k89olv$rr$1@newscl01ah.mathworks.com>...
> > > "Carl S." wrote in message <k885d7$h85$1@newscl01ah.mathworks.com>...
> > > > "Carl S." wrote in message <k879vh$k4n$1@newscl01ah.mathworks.com>...
> > > > > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h5rf05$nb8$1@fred.mathworks.com>...
> > > > > > Doug, I have a hard time to explain why using pointwise Neumann condition is dangerous in a simple term. So I continue to provide indication, and hope one of them will speak to you. If you are familiarize with boundary single layer potential, you will see that the normal derivative have a jump relation related to the local value of the potential (see Colton & Kress book). This show that we can set pointwise Neumann bc at any values, though the Banach L^p norm (1<=p<infinity) of the "trace" can be anything else. There are an infinity harmonic functions that meet the *point-wise* Neumann condition.
> > > > > >
> > > > > > 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
> > > > >
> > > > > Hi, I think this is the correct place to ask my question that is about Sobolev. The codes are ;
> > > > >
> > > > > 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 ?
> > > >
> > > > I tried this
> > > > SobolevGradEdgeTerm = (1./(eye(256)-gradient(u).^2)).*edgeTerm;
> > > > but I have still the same problem :(
> > >
> > > (eye(256)-gradient(u).^2) becomes zero. Therefore,
> > > 1./(eye(256)-gradient(u).^2)) becomes inf. Therefore, the active contour disappears
> > >
> > > How to solve this problem ?
> >
> > If I assign a small value that is 1e-10.*randi(1,1) instead of zeros then the above codes causes oversegmentation :(
>
> I have tried this
> temp2 = 1./(eye(size(u))-gradient(u).^2);
> temp2(isequal(temp2, inf)) = 1e-10;
> but temp2 have still inf values, why ? please help me

I tried this;
      [row,col] = find(temp2==Inf);
      temp2(row,col)=1e-10;
but, there is oversegmentation problem :(((((((
ImageAnalyst, where are you ? Any idea, any sugesstion ?

Subject: measured boundary conditions with pde toolbox

From: Dmitry

Date: 2 Apr, 2014 18:52:08

Message: 14 of 18

"Doug" wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> OK, I've answered my own question, and the answer is a boundary M-file.
>
> 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!

Hi Doug,
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...

Subject: measured boundary conditions with pde toolbox

From: Bill Greene

Date: 2 Apr, 2014 19:33:09

Message: 15 of 18

"Dmitry" wrote in message <lhhm8o$j16$1@newscl01ah.mathworks.com>...
> "Doug" wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> > OK, I've answered my own question, and the answer is a boundary M-file.
> >
> > 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!
>
> Hi Doug,
> 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...

There appears to be something wrong with your pdebound function or
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

Subject: measured boundary conditions with pde toolbox

From: Dmitry

Date: 3 Apr, 2014 04:56:08

Message: 16 of 18

"Bill Greene" wrote in message <lhholl$pur$1@newscl01ah.mathworks.com>...
> "Dmitry" wrote in message <lhhm8o$j16$1@newscl01ah.mathworks.com>...
> > "Doug" wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> > > OK, I've answered my own question, and the answer is a boundary M-file.
> > >
> > > 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!
> >
> > Hi Doug,
> > 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...
>
> There appears to be something wrong with your pdebound function or
> 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

Oh, thanks, Bill, I'd passed 'b' the wrong way...

Subject: measured boundary conditions with pde toolbox

From: Dmitry

Date: 7 Apr, 2014 15:11:08

Message: 17 of 18

"Dmitry" wrote in message <lhipl8$f2q$1@newscl01ah.mathworks.com>...
> "Bill Greene" wrote in message <lhholl$pur$1@newscl01ah.mathworks.com>...
> > "Dmitry" wrote in message <lhhm8o$j16$1@newscl01ah.mathworks.com>...
> > > "Doug" wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> > > > OK, I've answered my own question, and the answer is a boundary M-file.
> > > >
> > > > 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!
> > >
> > > Hi Doug,
> > > 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...
> >
> > There appears to be something wrong with your pdebound function or
> > 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
>
> Oh, thanks, Bill, I'd passed 'b' the wrong way...

Bill, could you help again please...
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?

Subject: measured boundary conditions with pde toolbox

From: Dmitry

Date: 8 Apr, 2014 10:00:12

Message: 18 of 18

"Dmitry" wrote in message <lhuf6c$ei2$1@newscl01ah.mathworks.com>...
> "Dmitry" wrote in message <lhipl8$f2q$1@newscl01ah.mathworks.com>...
> > "Bill Greene" wrote in message <lhholl$pur$1@newscl01ah.mathworks.com>...
> > > "Dmitry" wrote in message <lhhm8o$j16$1@newscl01ah.mathworks.com>...
> > > > "Doug" wrote in message <h5psrj$8g9$1@fred.mathworks.com>...
> > > > > OK, I've answered my own question, and the answer is a boundary M-file.
> > > > >
> > > > > 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!
> > > >
> > > > Hi Doug,
> > > > 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...
> > >
> > > There appears to be something wrong with your pdebound function or
> > > 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
> >
> > Oh, thanks, Bill, I'd passed 'b' the wrong way...
>
> Bill, could you help again please...
> 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?

Thanks a lot!
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.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us