Thread Subject: Problem solving bvp of Meniscus

I've been trying to solve the following equation:
y''/(1+y'^2)^(3/2) = y + a

and a is a constant value that should be determined after solving the equation.

That's a boundary value problem, with the conditions:

z(infinity) = 0, z'(infinity) = 0 and z'(0) = -cotg(theta)

I wrote the equation into matlab and a file named ode_men.m with:

function dydx = ode_men(x, y, a)
    % y(1) = y
    dydx = [y(2) % = y'
            (y(1)+a)*(1+y(2)^2)^(3/2)]; % = y''

The initial guess was an exponential guess, a good guess I think, with:

function y = guess(x)
    a = 0.39;
    h = a*sqrt(1-sin(pi/6));
    b = -1*(log(10^(-6)/h))/1000;
    y = [h*exp(-b*x)
         -b*h*exp(-b*x)];
 in a file guess.m

My boundary values are in a file named bvp.m :

function res = bvp(ya, yb, a)
    theta = (pi/180)*30; %thirty degrees
    res = [ya(2) % z(infinity)
           yb(2) % z'(infinity)
           yb(1) + cot(theta)]; % initial angle cotg(theta)

to solve, I used:

solinit = bvpinit(linspace(0,1,10), @guess, 1);
% interval in x, values in y for each x and then the guess for a
sol = bvp4c(@ode_men, @bvp, solinit)

my return was:

sol =

        solver: 'bvp4c'
             x: [0 0.0556 0.1111 0.1667 0.2222 0.2778 0.3333 0.3889 0.4444 0.5000 0.5556 0.6111 0.6667 0.7222 0.7778 0.8333 0.8889 0.9444 1]
             y: [2x19 double]
            yp: [2x19 double]
    parameters: 1.7321
         stats: [1x1 struct]

The problem, it's returning a constant function. Totally different of the expected result.

>> sol.y

ans =

  Columns 1 through 16

   -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321 -1.7321
         0 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000

  Columns 17 through 19

   -1.7321 -1.7321 -1.7321
    0.0000 -0.0000 -0.0000

These values -1.7321 are the values of y for each x. And the zeros are values of y' for each x.
My a value returned was 1.7321, and it's equal to all the values of y. I think it came from the third line of my bvp that I wrote yb(1) + cot(theta) corresponding to y'(0) = -cotg(30°) = -sqrt(3).

Changing to boundary values and making doing a similar guess, I had other errors, some with jacobian. I would like to solve this equation for that boundary conditions I already said and for another:

y(0) = 0, y'(0) = 0 and y'(1) = cotg(theta)

This last one, I was able to solve with maple. But the first one, that I wrote the conditions before, didn't work.


Can somebody help me? Am I doing something wrong?

Thanks for all,
Regards

Ataías Pereira Reis <ataiasreis@gmail.com> wrote in message <1207681.241.1319819354868.JavaMail.geo-discussion-forums@yqgn17>...
> I've been trying to solve the following equation:
> y''/(1+y'^2)^(3/2) = y + a
>
> and a is a constant value that should be determined after solving the equation.
>
> That's a boundary value problem, with the conditions:
>
> z(infinity) = 0, z'(infinity) = 0 and z'(0) = -cotg(theta)
> .........
 - - - - - - - - - -
  If by z(infinity) = 0, z'(infinity) = 0 you mean that these conditions apply to y = z, then you would also need y"(infinity) = 0. Any other value would make y'(infinity) = 0 impossible. This in turn implies that the quantity 'a' must be zero to satisfy your original differential equation. Hence you have just the equation

  y" = y*(1+y'^2)^(3/2)

  This is not really a boundary value problem if you do the following. Define

 v = y' = dy/dt

Then

 y" = dv/dt = (dv/dy)*(dy/dt) = dv/dy * v = 1/2*dv^2/dy

which gives you the equation

 1/2*dv^2/dy = y*(1+v^2)^(3/2)

or written as differentials

 1/2*(1+v^2)^(-3/2)*dv^2 = y*dy

Integrating this gives

 -(1+v^2)^(-1/2) = y^2/2 + C

where C is the constant of integration. The condition y(infinity) = 0 and v(infinity) = 0 then implies that this constant must be C = -1, giving

 y'^2 = v^2 = 1/(1-y^2/2)^2 - 1 = (y^2-y^4/4)/(1-y^2/2)^2

or

 y' = (+ or -) sqrt((y^2-y^4/4)/(1-y^2/2)^2)

This is then the equation you hand to one of matlab's ode functions.

  The only additional condition you need now is y(0) or y'(0) and the sign will immediately determine the sign needed in the equation for y'.

Roger Stafford

"Roger Stafford" wrote in message <j8equl$em1$1@newscl01ah.mathworks.com>...
> .......
> y' = (+ or -) sqrt((y^2-y^4/4)/(1-y^2/2)^2)
> .......
- - - - - - -
  An afterthought. If you write that last equation like this, there will be no uncertainly as to the sign:

 y' = y/(y^2-2)*sqrt(4-y^2)

  Note that you must always have -sqrt(2) < y < +sqrt(2), and y' and y must always be of opposite signs, (given the assumed asymptote at infinity.)

Roger Stafford

Hello Roger! First off, thank you for your effort in your answer. Anyway, it's not supposed that a=0, because it's related with the height of the meniscus. By the way, I think you didn't considered the fact of y'(0) = -cotg(theta).

If a should be zero in this equation, it's a wrong equation. I was trying to solve another equation before, but I have a lot of errors, the graph had never fitted the expected. This equation was:

\frac{2y}{a^2} - \frac{z''}{(1+(z')^2)^{3/2}} = 0

I obtained some graphs with this equation (terrible graphics), using the conditions:

function res = bvp(ya, yb, a)
    theta = (pi/180)*30; %thirty degrees
    res = [ya(2) % z(infinity)
           yb(2) % z'(infinity)
           yb(1) + cot(theta)]; % initial angle cotg(theta)
that I'd shown before.

Well, in this equation, a can not be zero, as it's shown as a divisor. The problem with the graphs is that I obtained a lot of horizontal lines (not what I wanted) or an exponential with a very high value when x = 0. Not the case, because the height is the value of y(0), and it's not very big. The a value should be something greater than 0 and less than 1. The height is obtained knowing the a value, but I don't remember the equation for h.

Well, does anyone have an idea how to fix the bug of a horizontal line? I don't know why this happens (well, I showed the returned values in the first post, still that problem).

Thank you for all,
Regards

Ataías Pereira Reis <ataiasreis@gmail.com> wrote in message <24161917.1073.1320078271146.JavaMail.geo-discussion-forums@yqp37>...
> Hello Roger! First off, thank you for your effort in your answer. Anyway, it's not supposed that a=0, because it's related with the height of the meniscus. By the way, I think you didn't considered the fact of y'(0) = -cotg(theta).
>
> If a should be zero in this equation, it's a wrong equation. ........
- - - - - - - - -
  The original differential equation you presented was:

 y"/(1+y'^2)^(3/2) = y + a

However you stated conditions in terms of a variable z. I assumed previously that you meant that y and z were to be regarded as the same variable. If so, then your requirement that z(infinity) = 0 in the case of your equation will force all derivatives at infinity to also be zero and this does imply that 'a' must be zero. Just plug 0 into the above equation for y, y', and y" and see what you get. What we are talking about here is a requirement that y approach zero asymptotically which is far different from merely requiring a particular value of y at a finite point. It is a strong requirement.

  Please state whether z and y are understood to be the same quantities. If not, how do you define z? It would be helpful if you could restate your conditions on y unambiguously without a reference to z.

  In any case, since the underlying independent variable t is not explicitly present in your equation, the trick I showed you of replacing y" by 1/2*dv^2/dy where v = y' can still be used in your equation regardless of what your boundary conditions are and what the value of 'a' is. It reduces your equation to first order and that may be useful in solving it.

Roger Stafford