Solving Second-Order BVP with unknown Constants using fsolve

7 Ansichten (letzte 30 Tage)
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 Kommentare
Torsten
Torsten am 7 Okt. 2024
Use "bvp4c". It can solve boundary value problems with unknown parameters.
Mark Fasano
Mark Fasano am 7 Okt. 2024
Bearbeitet: Mark Fasano am 7 Okt. 2024
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)

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
Torsten am 7 Okt. 2024
Bearbeitet: Torsten am 8 Okt. 2024
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 Kommentare
Mark Fasano
Mark Fasano am 7 Okt. 2024
Bearbeitet: Mark Fasano am 7 Okt. 2024
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
Torsten am 8 Okt. 2024
Bearbeitet: Torsten am 8 Okt. 2024
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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte


Version

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by