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

# Thread Subject: Problem solving bvp of Meniscus

 Subject: Problem solving bvp of Meniscus From: AtaÃ­as Pereira Reis Date: 28 Oct, 2011 16:29:14 Message: 1 of 5 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
 Subject: Problem solving bvp of Meniscus From: Roger Stafford Date: 28 Oct, 2011 18:06:13 Message: 2 of 5 AtaÃ­as Pereira Reis 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
 Subject: Problem solving bvp of Meniscus From: Roger Stafford Date: 29 Oct, 2011 00:46:32 Message: 3 of 5 "Roger Stafford" wrote in message ... > ....... > 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
 Subject: Problem solving bvp of Meniscus From: AtaÃ­as Pereira Reis Date: 31 Oct, 2011 16:24:31 Message: 4 of 5 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