Solving Second-Order BVP with unknown Constants using fsolve

I am trying to solve the following problem numerically using Matlab's fsolve function. The problem is specified by:

with the boundary conditions that

, ,

The unknowns in the problem are the height profile , the angle , and the unknown constant p. The problem is in fact well posed because we have the two "extra" boundary conditions are required for us to solve for the two unknown constants . I devised the following numerical scheme:

[Eq. **]

where and the length of the domain is . The number of points used to discretize space is N. The above scheme is to be applied at . At , we have to be mindful of the boundary conditions. I used the contact angle BCs to come up with equations for the "ghost points", or points just outside the domain. These are given by:

, and

Utilizing these ghost point equations and the fact that , the equations for are given by:

,

My idea is to now solve Eq. [**] for with the four equations immediately above. This gives a total of equations for the unknowns which are: . To implement this in Matlab, I wanted to use an iterative method which calls fsolve at each iteration. I, however, got a bit confused as to the best way of defining these equations as function handles to then use in the fsolve command. Below is my definition of the function that should define the equations:

function F = finDiffScheme(hold,N) % here h is the value of h at the previous iteration level
    xf = 6.17415;
    L = 2*xf; dx=L/N; x=-xf:dx:xf;
    theta = 39*(pi/180);
    S = 5; ki=.37; alpha1 = 2.69;
    F{1} = @(y) 2*y(1) - 2*dx*tan(theta-y(N+1))-dx^2*S+dx^2*y(N);
    F{2} = @(y) y(2)-(2+dx^2)*y(1)-dx^2*S*exp(-2*ki*((x(2)+xf)+alpha1*hold(1)))+dx^2*y(N);
    for i = 2:N-2
        F{i+1} = @(y) y(i+1)-(2+dx^2)*y(i)+y(i-1)-dx^2*S*exp(-2*ki*((x(i+1)+xf)+alpha1*hold(i)))+dx^2*y(N);
    end
    F{N} = @(y) -(2+dx^2)*y(N-1)+y(N-2)-dx^2*S*exp(-2*ki*((x(N)+xf)+alpha1*hold(N-1)))+dx^2*y(N);
    F{N+1} = @(y) 2*y(N-1) - 2*dx*tan(theta+y(N+1))-dx^2*S*exp(-4*ki*xf)+dx^2*y(N);
end

If I run this code, then the output is a cell array called F which is size 1 x N+1. My concern here is that if I look at F, it seems that the numerical values of the known parameters S, ki, alpha1, dx, xf are all ignored. This is very confusing to me since F should only be a function of y not also those parameters. Is there a reason for this? For reference, I call this function in the script below:

clear; clc; 
% Set parameters first: xf, theta
xf = 6.17415;
theta = 39*(pi/180);
h0 = 0;
hN = 0;
iters = 3;
N=256; L=2*xf; dx = L/N; x=-xf:dx:xf;

uold = zeros(1,N+3);   % First n-1 elements are h_1, ..., h_{n-1}; Last 2 elements are p and \theta_d
hold1 = (1 + (sinh(x-xf)-sinh(x+xf))/sinh(2*xf))*.809827;
uold(1,1:N+1) = hold1;
uold(1,N+2:N+3) = [.809827,0];
plot(x,hold1)
hold on
F = finDiffScheme(hold1,N);
F2 = @(z) cellfun(@(y) y(z), F)

for j=1:iters
    u = fsolve(F2, uold)
    plot(x,u(1,1:N+1))
    disp(u(1,N+2:N+3))
    uold = u;
end

I am a bit concerned that the fsolve command is somehow also solving for these parameters. What exactly am I doing wrong here?

For reference, I was trying to define my solution vector as so that for , , and .

[Torsten]

Use "bvp4c". It can solve boundary value problems with unknown parameters.

https://uk.mathworks.com/help/matlab/ref/bvp4c.html

[Mark Fasano]

Thanks for pointing me in that direction. I am trying to implement this as directed by the example in this link where they use bvp4c to solve a second order system with one unknown parameter. I am, however, getting errors that confuse me. The error message reads: "Not enough input arguments. Error in solution>mat4bc (line 22): ya(2)-tan(.6807-thetad)." I honestly have no idea what to make of this. I thought I followed that example exactly however am a bit confused. My code is below.

xmesh = linspace(-6.17415,6.17415,64);
solinit = bvpinit(xmesh,@mat4init,[.809827,0]);
sol = bvp4c(@mat4ode,@mat4bc, solinit);

xint = linspace(-6.17415,6.17415);
Sxint = deval(sol,xint);

plot(xint,Sxint)
axis([-6.17415 6.17415 -4 4])
xlabel('x')
ylabel('h')
legend('h','h''')

function dydx = mat4ode(x,y,p,thetad) % Equation being solved
    dydx = [y(2)
        y(1)+2*exp(-2*.37*((x+6.17415)+2.69*y(1)))-p];
end

function res = mat4bc(ya,yb,p,thetad) % Boundary conditions
    res = [ya(1) 
           yb(1)
           ya(2) - tan(.6807-thetad)
           yb(2) + tan(.6807+thetad)];
end

function yinit = mat4init(x) % Initial guess
    yinit = [(1 + (sinh(x-6.17415)-sinh(x+6.17415))/sinh(2*6.17415))*.809827
        .809827*(cosh(x-6.17415)-cosh(x+6.17415))/sinh(2*6.17415)];
end

I started by defining and so that the ODE can be written as the first-order system:

I then coded this up into the function mat4ode. Notice that I manually inserted the numerical value for the known parameters in my function - that is, I used .

I then rewrote my boundary conditions in terms of . The Dirichlet BCs are written first in the function mat4bc and then the Neumann conditions last. I was a bit confused as to what exactly ya and yb are. I believe, in my case and in the example, they are vectors of size 2 that contain the values in the vector ya and similarly for the vector yb (please let me know if this is incorrect).

Finally I defined my initial guess. I believe this is supposed to be a vector of sized 2 which holds an initial guess for and in that order. In my case, I just used the exact solution to the problem when .

It seems that the issue is with the boundary conditions given the error message. I wonder if there is a more fundamental issue given the fact that doesn't actually appear in the ODE but only the BCs (similarly that p doesn't appear in the BC but appears in the ODE)

Answer 1

[Torsten]

solinit = bvpinit(linspace(-6.17415,6.17415,64),[0;0],[.809827,0]);
sol = bvp4c(@mat4ode,@mat4bc, solinit);

xint = linspace(-6.17415,6.17415);
Sxint = deval(sol,xint);

A = trapz(xint,Sxint(1,:))

A = 12.1055

plot(xint,Sxint)
axis([-6.17415 6.17415 -4 4])
xlabel('x')
ylabel('h')
legend('h','h''')

function dydx = mat4ode(x,y,p)
    dydx = [y(2)
        y(1)+2*exp(-2*.37*((x+6.17415)+2.69*y(1)))-p(1)];
end

function res = mat4bc(ya,yb,p)
    res = [ya(1)
           yb(1)
           ya(2) - tan(.6807-p(2))
           yb(2) + tan(.6807+p(2))];
end

[Mark Fasano]

Thank you for the help! I notice two differences between my work and your solution:

1. It seems you combined the unknown parameters into a vector p. Is this always necessary? i.e. imagine a problem with three unknown parameters would you just make p a vector of sized 3?
2. I also noticed you called bvpinit differently. What exactly was wrong with the way I initially did it?

[Torsten]

It seems you combined the unknown parameters into a vector p. Is this always necessary? i.e. imagine a problem with three unknown parameters would you just make p a vector of sized 3?

Yes.

I also noticed you called bvpinit differently. What exactly was wrong with the way I initially did it?

Nothing, I just didn't find your init-function. You must have added it later. The code will also work with your complicated initial guesses for h and h' and will give the same result.

[Mark Fasano]

Is there anyway to incorporate integral constraints into this system? One thing I want to try to incorporate as well is the constraint that we know the area of the drop, i.e. that for some known constant A we have:

I saw this post which described one way to indirectly do this. To implement this, I tried to add the extra equation

and then to also include the boundary conditions

To do this, I took the code in the answer above and just changed the following functions:

function dydx = mat4ode(x,y,p)
    dydx = [y(2)
        y(1)+2*exp(-2*.37*((x+6.17415)+2.69*y(1)))-p(1)
        y(1)];
end
function res = mat4bc(ya,yb,p)
    res = [ya(1)
           yb(1)
           ya(2) - tan(.6807-p(2))
           yb(2) + tan(.6807+p(2))
           ya(3) - 8
           yb(3)];
end

and then initially just ensuring the initial guess for y is a vector of length 3. The problem with this is I get the following error: "The boundary condition function BCFUN should return a column vector of length 5." This intuitively makes sense to me as I am trying to solve a system of 3 first-order ODEs with 2 unknown parameters so I should only have 5 degrees of freedom. Is there no way around this issue.

Since I can't figure out how to incorporate that into the code, I alternatively thought I could have the script solve for (the endpoint of the domain) given some value of A (say 8 as in the example above). This, however, doesn't seem possible using this function considering you wouldn't even be able to create your xmesh. Is that correct? Do you know what function I should use in that case?

[Torsten]

The quantity

A = integral_{x=-xf}^{x=xf} h(x) dx 

is a deduced quantity from the solution for h, i.e. it is fixed once you have solved your original problem. You can't prescribe a value for A.

If you want to prescribe a value for A, you have to give up some condition of your original problem (maybe the value of xf ? Or some boundary condition ? ).

I added the computation of the deduced quantity A to the code above once the solution for h is available.