Solving Second-Order BVP with unknown Constants using fsolve
7 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Mark Fasano
am 7 Okt. 2024
Bearbeitet: Torsten
am 8 Okt. 2024
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
Akzeptierte Antwort
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,:))
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
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.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!