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 . 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 . 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 .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 . 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 . 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 . 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 . 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 . 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 . ode, fsolve, functions, cellfun MATLAB Answers — New Questions
