Solving Second-Order BVP with unknown Constants using fsolve

Question

Mark Fasano 2024-10-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2157990-solving-second-order-bvp-with-unknown-constants-using-fsolve

编辑： Torsten 2024-10-8

采纳的回答： Torsten

在 MATLAB Online 中打开

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

.

2 个评论
显示无隐藏无

Torsten 2024-10-7

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

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

Mark Fasano 2024-10-7

编辑：Mark Fasano 2024-10-7

在 MATLAB Online 中打开

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 2024-10-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2157990-solving-second-order-bvp-with-unknown-constants-using-fsolve#answer_1528250

编辑：Torsten 2024-10-8

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mark Fasano 2024-10-7

编辑：Mark Fasano 2024-10-7

在 MATLAB Online 中打开

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 2024-10-8

编辑：Torsten 2024-10-8

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Solving Second-Order BVP with unknown Constants using fsolve

2 个评论
显示无隐藏无

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving Second-Order BVP with unknown Constants using fsolve

2 个评论 显示 无隐藏 无

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

4 个评论
显示 2更早的评论隐藏 2更早的评论