Solving a 6th order non-linear differential equation in Matlab
Show older comments
Hello everyone,
I am trying to solve a high-order non linear differential equation presented right below with 
:
:
As boundary conditions, I have the following ones :
I am trying to apply a shooting method algorithm in order to solve this equation on 
, I have converted it in a system of 1-st order differential equations as following :
, I have converted it in a system of 1-st order differential equations as following :
with : 
According to the boundary conditions we thus have : 
and we have to guess the remaining ones : 
and we have to guess the remaining ones : I have attached to my post, the algorithms that I wrote which are inspired by biblographical researchs : the main code is called "shoot4nl.m" and I have tried to guess some initial values in order to run the algorithm. Also, I have choses a values of 
very "high" in order to match the conditions...
very "high" in order to match the conditions...By doing so, I obtain the following error :
Warning: Matrix is singular to working precision.
> In newtonRaphson2 (line 23)
In shoot4nl (line 14)
Error using newtonRaphson2
Too many iterations
Error in shoot4nl (line 14)
u = newtonRaphson2(@residual,u);
I assume that means that there is a division by 0 that occurs in the process but I remain without any idea to solve this problem...
Could someone help me or bring his mathematical expertise in order to show me some hints ?
Thank you in advance,
Best regards.
2 Comments
Bruno Luong
on 7 Apr 2022
Edited: Bruno Luong
on 7 Apr 2022
6th order???? That's quite large and you need at least to be careful how is the scaling of your function, for axample if you have eta that is large/small compare to unit (1), the order of state vector f' has exponetially large discrepency, and the Jacobian of the syetem becomes numerically ill-condition.
Try to reformulate ODE of
f(eta)
as ODE of
g(xi) := f(xi*etascale)
where etascale is a typical scale of eta. You might get better chance to solve your system.
Wissem
on 7 Apr 2022
Answers (2)
I think that nobody in this forum will dive into selfmade software if there are well-tested MATLAB alternatives.
So these few lines of code can get you started - backed with professional software:
bc0 = [1 1 1 1];
bc = fsolve(@fun,bc0);
F0 = [1e-5;0;bc(1);bc(2);bc(3);bc(4)]
etaspan = [-10,10];
[eta,F] = ode45(@fun_ode,etaspan,F0);
plot(eta,F(:,1))
function res = fun(bc)
F0 = [1e-5;0;bc(1);bc(2);bc(3);bc(4)];
etaspan = [-10, 0, 10];
[eta,F] = ode45(@fun_ode,etaspan,F0);
res(1) = F(2,2);
res(2) = F(2,4);
res(3) = F(3,1);
res(4) = F(3,2);
end
function dFdeta = fun_ode(eta,F)
dFdeta = [F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
30 Comments
Wissem
on 6 Apr 2022
Torsten
on 6 Apr 2022
res defines the residuals of the shooting method at eta = 0 and eta = eta0:
res(1) = F(2,2) is the deviation of F' from 0 at eta = 0
res(2) = F(2,4) is the deviation of F''' from 0 at eta = 0
res(3) = F(3,1) is the deviation of F from 0 at eta = eta0
res(4) = F(3,2) is the deviation of F' from 0 at eta = eta0
fsolve tries to adjust the initial conditions for F'', F''', F'''', F''''' at eta = -eta0 such that prescribed conditions at eta = 0 and eta = eta0 hold.
So different from your setting, integration starts at eta = -eta0, not at eta = 0.
Wissem
on 6 Apr 2022
Torsten
on 6 Apr 2022
I used octave and don't get an error.
Why do you want to adjust stepsize ? The solver usually selects the adequate stepsize in order to satisfy the prescribed error tolerance. So if you lower the error tolerance, stepsize will be larger and vice versa.
Wissem
on 6 Apr 2022
Torsten
on 6 Apr 2022
Adjusting the stepsize won't help here.
The crucial adjustments are the etaspan over which you want to integrate and the value for F that you prescribe at eta = eta0.
A value F=0 is not possible since you must divide by F.
The following code works,e.g., and gives the included solution for F:
bc0 = [2.2504e+06 -2.5004e+05 -2.2123e+06 5.6616e+05]
etaspan = [-30 30];
bc = fsolve(@(bc)fun(bc,etaspan),bc0);
bc_total = [1e-7;0;bc(1);bc(2);bc(3);bc(4)];
options = odeset('InitialStep',1e-9);%,'MaxOrder',1);
[eta,F] = ode45(@fun_ode,etaspan,bc_total);%,options);
plot(eta,F(:,1))
function res = fun(bc,etaspan)
bc_total = [1e-7;0;bc(1);bc(2);bc(3);bc(4)];
etaspan_ode = [etaspan(1),0,etaspan(2)];
options = odeset('InitialStep',1e-9);%,'MaxOrder',1);
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);%,options);
res(1) = F(2,2);
res(2) = F(2,4);
res(3) = F(3,1)-1e-7;
res(4) = F(3,2);
end
function dFdeta = fun_ode(eta,F)
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
Wissem
on 7 Apr 2022
Could you please take the time to explain the modifications that you provided from your previous correction ?
I just made the problem harder to solve: F(-eta0) = F(eta0) = 1e-7 instead of 1e-5 and the value for eta0 is bigger now.
I guess the boundary condition at eta0 is given at +/-Inf ?
Maybe mapping (-oo,oo) to [-1:1] and integrating over [-1:1] could make the problem easier to solve. A possible coordinate transformation is eta~ = atanh(eta).
Also using the finite-difference method instead of shooting (i.e.using bvp4c) could make the solution process more stable.
Wissem
on 7 Apr 2022
Torsten
on 7 Apr 2022
If you use my code under octave, it should work. I don't think that something has been changed in the fsolve or ode45 functions. Which version do you use ? My version is 6.4.0.
I tried to set up the problem for bvp4c, but I couldn't test it - so there might be erros as well as that the solver doesn't succeed.
The code should somehow look like
eta0 = 30;
etamesh = [linspace(-eta0,0,100),linspace(0,eta0,100)];
Finit = [1.0e-7; 0; 0; 0; 0; 0];
solinit = bvpinit(etamesh,Finit)
sol = bvp4c(@fun_ode, @bc, solinit);
plot(sol.x,sol.y(1,:))
function dFdeta = fun_ode(eta,F,region)
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
function res = bc(FL,FR)
res = [FL(1,1) - 1.0e-7 ; FL(2,1);...
FR(2,1) ; FR(4,1);...
FR(1,1)-FL(1,2) ; FR(2,1)-FL(2,2) ; FR(3,1)-FL(3,2) ; FR(4,1)-FL(4,2) ; FR(5,1)-FL(5,2) ; FR(6,1)-FL(6,2);...
FR(1,2) - 1.0e-7 ; FR(2,2)];
end
Wissem
on 7 Apr 2022
Torsten
on 7 Apr 2022
Use the code below, name it "main.m", load it in the octave editor and run it.
function main
bc0 = [2.2504e+06 -2.5004e+05 -2.2123e+06 5.6616e+05]
etaspan = [-30 30];
bc = fsolve(@(bc)fun(bc,etaspan),bc0);
bc_total = [1e-7;0;bc(1);bc(2);bc(3);bc(4)];
options = odeset('InitialStep',1e-9);%,'MaxOrder',1);
[eta,F] = ode45(@fun_ode,etaspan,bc_total);%,options);
plot(eta,F(:,1))
end
function res = fun(bc,etaspan)
bc_total = [1e-7;0;bc(1);bc(2);bc(3);bc(4)];
etaspan_ode = [etaspan(1),0,etaspan(2)];
options = odeset('InitialStep',1e-9);%,'MaxOrder',1);
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);%,options);
res(1) = F(2,2);
res(2) = F(2,4);
res(3) = F(3,1)-1e-7;
res(4) = F(3,2);
end
function dFdeta = fun_ode(eta,F)
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
Wissem
on 8 Apr 2022
Why did you choose to compute 4 residuals in the function ?
As I already tried to explain to you, I start the ode-integrator with the two conditions set at eta = -eta0 (F = F' = 0) and estimates for the remaining four initial conditions for F'',F''',F'''' and F'''''. You defined four other conditions (two at 0, two at eta = eta0 for F) that should hold. So we have four values that we have to nullify, namely F' = F''' = 0 at eta=0 and F = F' = 0 at eta=eta0. These are the four equations that fsolve have to solve. In the res vector, the errors are returned to fsolve and it is prompted to supply better guesses for F''(-eta0),...,F'''''(-eta0) in order to nullify these errors.
- I understand the aim of the line :
bc = fsolve(@(bc)fun(bc,etaspan),bc0);
However, I am not familiar with this syntax (anonymous function) : why did you choose to write it this way ?
"etaspan" is needed in "fun" to define the integration interval for the ode integrator. The line
bc = fsolve(@(bc)fun(bc,etaspan),bc0);
has the effect that "etaspan" is passed to "fun" as an additional parameter that can be used in the function. If you don't like this, you will have to change "etaspan" both in the calling function and in "fun" if you make changes to eta0.
I have strong hope that "bvp4c" can solve your problem more stable than the ad-hoc shooting method I supplied. Did you already test the code I gave you in MATLAB ?
Wissem
on 8 Apr 2022
Wissem
on 8 Apr 2022
@Torsten Do you have any idea about the fact that the same code runs well on octave but keeps telling me there is an error on fsolve in MATLAB ??
No solution found.
fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.
Wissem
on 8 Apr 2022
Since I don't have MATLAB available, I can't test.
But this setting (eta0 = 30, F(eta0)=1e-7) is already very difficult for the solvers. Start with less aggressive values (eta0 = 10, F(eta0) = 1e-3 or something). If you succeed, use the solution as initial value for making eta0 larger and F(eta0) smaller.
But also under OCTAVE, the calculations were quite unstable and small changes in eta0 and F(eta0) led to big changes in the height at eta=0. Although you got the error message of a singular Jacobian, you should try bvp4c or bvp5c with the code I gave to you. Solving boundary value problems with finite differencing is in general much more stable than shooting.
Also solving on a finite interval (e.g. [-1:1] by the coordinate transformation eta~ = tanh(eta)) could stabilize the computations.
Somehow I have the feeling that a condition is missing in the problem formulation. If you start integration at x=0, you have three conditions (because I think the solution is symmetric to eta=0), namely F'(0)=F'''(0)=F'''''(0)=0. Additionally, you have the two conditions at Infinity: F(infinity) = F'(Infinity) = 0. But what is the necessary 6th condition ?
Do you think it's the good numerical reformulation for this particular script ?
Are you aware that this is a totally different problem ? The results of the old formulation gave a value in the order of 1e6 in eta = 0. Now you prescribe F = 1 in eta = 0. Do you think that both formulations have anything in common except for the underlying differential equation ?
This should be the correct code for the reformulated problem, but it doesn't converge:
bc0 = [-2.5004e+05 -2.2123e+06 5.6616e+05]
etaspan = [0 10];
bc = fsolve(@(bc)fun(bc,etaspan),bc0); % Adjusting the initial conditions for F''',F'''',F''''' at eta = eta0
bc_total = [1;0;bc(1);0;bc(2);bc(3)];
[eta,F] = ode45(@fun_ode,etaspan,bc_total);
plot(eta,F(:,1))
% Creating residuals for the shooting method
function res = fun(bc,etaspan)
bc_total = [1;0;bc(1);0;bc(2);bc(3)];
etaspan_ode = [etaspan(1),etaspan(2)];
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);
res(1) = F(end,1)-1e-2; % Deviation of F from 0 at eta = eta0
res(2) = F(end,2); % Deviation of F' from 0 at eta = eta0
res(3) = F(end,3); % Deviation of F'' from 0 at eta = eta0
end
% System of 1-st order ODE
function dFdeta = fun_ode(eta,F)
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
Wissem
on 9 Apr 2022
@Torsten ,
I have made a mistake in my previous post : I am following your advice trying to solve the problem on the 
interval only. In this case, I have introduced a new boundary condition such as :
interval only. In this case, I have introduced a new boundary condition such as : 
I think that this is the only way to solve the problem because as you noticed : the solution obtained with OCTAVE was really unstable. The solution value at eta = 0 is very important so I need a robust way to determine it.
Here are the modifications provided to the script :
% Main script for the shooting method
bc0 = [2.2504e+06 2.2123e-6 5.6616e-05];
etaspan = [0 10];
bc = fsolve(@(bc)fun(bc,etaspan),bc0); % Adjusting the initial conditions for F''',F'''',F''''' at eta = eta0
bc_total = [bc(1);0;bc(2);0;bc(3);0];
[eta,F] = ode45(@fun_ode,etaspan,bc_total);
plot(eta,F(:,1))
% Creating residuals for the shooting method
function res = fun(bc,etaspan)
bc_total = [bc(1);0;bc(2);0;bc(3);0];
etaspan_ode = [etaspan(1),etaspan(2)];
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);
res(1) = F(end,1)-1e-2; % Deviation of F from 0 at eta = eta0
res(2) = F(end,2); % Deviation of F' from 0 at eta = eta0
res(3) = F(end,3); % Deviation of F'' from 0 at eta = eta0
end
% System of 1-st order ODE
function dFdeta = fun_ode(eta,F)
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
With always an error...
I will also try my old formulations (on 
) but with 
: this is a quiet difficult problem to solve...
) but with : this is a quiet difficult problem to solve...Thank you again.
The solution value at eta = 0 is very important so I need a robust way to determine it.
I don't want to discourage you, but experimenting with the old code showed that the value at eta=0 is strongly dependent on the eta0 chosen.
This works under OCTAVE:
bc0=[ 3.8309 0.4494 -0.1066];
etaspan = [0 12];
bc = fsolve(@(bc)fun(bc,etaspan),bc0); % Adjusting the initial conditions for F''',F'''',F''''' at eta = eta0
bc_total = [bc(1);0;bc(2);0;bc(3);0];
options = odeset('InitialStep',1e-9);
[eta,F] = ode45(@fun_ode,etaspan,bc_total,options);
bc
plot(eta,F(:,1))
end
% Creating residuals for the shooting method
function res = fun(bc,etaspan)
bc_total = [bc(1);0;bc(2);0;bc(3);0];
etaspan_ode = [etaspan(1),etaspan(2)];
options = odeset('InitialStep',1e-9);
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total,options);
res(1) = F(end,1)-1e-5; % Deviation of F from 0 at eta = eta0
res(2) = F(end,2); % Deviation of F' from 0 at eta = eta0
res(3) = F(end,3); % Deviation of F'' from 0 at eta = eta0
end
% System of 1-st order ODE
function dFdeta = fun_ode(eta,F)
%F
dFdeta=[F(2);F(3);F(4);F(5);F(6);(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
Torsten
on 9 Apr 2022
bvp4c or bvp5c is not part of the OCTAVE software.
This examples shows how to set up problems where intermediate conditions have to be set at inner grid points (like at eta=0 for the integration over the interval [-eta0,eta0]):
"bvpinit" is no package - it's just a function that is called to set initial conditions for the BVP. It doesn't help you to do so - it just offers you the possibility to do it as a constant value or via a function where you can set the values dependent on eta.
Wissem
on 9 Apr 2022
I suggest you start with the solution obtained from the shooting method as initial guess.
I'm curious whether bvp4c will reproduce the solution.
Using bvp4c
eta0 = 10;
eta = linspace(0,eta0,100);
yinit = [1e4;0;0;0;0;0];
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
plot(eta,F);
function dF = fun_ode(eta,F)
dF=[F(2);
F(3);
F(4);
F(5);
F(6);
(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
function res = bcfun(ya,yb)
res =[ya(2); % F'(0)
ya(4); % F'''(0)
ya(6); % F'''''(0)
yb(1:3)]; % F(eta0), F'(eta0), F''(eta0),
end
79 Comments
The value of F at eta0 is 625, not 0.
Most probably,
res =[ya(2); ya(4); ya(6); yb(1); yb(2); yb(3)];
instead of
res =[ya(2); % F'(0)
ya(4); % F'''(0)
ya(6); % F'''''(0)
yb(1:3)]; % F(eta0), F'(eta0), F''(eta0),
@Bruno Luong @Torsten : Thank you again for your contributions. However, I am a bit confused : one of the boundary condition imposed states that 
but here the plotting solution gives 625 which is not in adequation with the analytical formulation ?? The solution must surely be a bell shape centered on 0 as obtained in @Torsten previous plots.
but here the plotting solution gives 625 which is not in adequation with the analytical formulation ?? The solution must surely be a bell shape centered on 0 as obtained in @Torsten previous plots.Thank you again,
Torsten
on 10 Apr 2022
Is the result the same if you run Bruno's code with
res =[ya(2); % F'(0)
ya(4); % F'''(0)
ya(6); % F'''''(0)
yb(1:3)]; % F(eta0), F'(eta0), F''(eta0),
replaced by
res =[ya(2); ya(4); ya(6); yb(1); yb(2); yb(3)];
?
Bruno Luong
on 10 Apr 2022
Edited: Bruno Luong
on 10 Apr 2022
I think it converges to a local minima depending on yinit.
I don't know the math/stability behind those equations, but if you replace with Lorenz equation, he showed that the bvp is impossible to solve numerically since it is chaotic the dependency of the final state to the initial state.
Wissem
on 10 Apr 2022
@Torsten Yes, the result remains the same : we obtain the same plot with 
instead of 0 and I don't understand why since we are imposing it as a boundary condition.
instead of 0 and I don't understand why since we are imposing it as a boundary condition.@Bruno Luong Thank you for your comment : but how did you came up with this conclusion because, I think that the bvp solver uses Lobatto integration method isn't it ?
Thank you.
Bruno Luong
on 10 Apr 2022
Edited: Bruno Luong
on 10 Apr 2022
It's not problem off Lobatto or integration, if your ode equation is non linear chance that the objective function (square norm of res) is non convex and it likely stuck at the local minimum.
Solving bcp non-linear 6th order ode is insane, unless you know inside-out the physics and maths of the system you are trying to solve, I wouldn't expect it works with 10 lines of careless code (mine).
Bruno Luong
on 10 Apr 2022
Edited: Bruno Luong
on 10 Apr 2022
No the ode seems to be alright, the system doesn't seem to be stiff, if you solve Cauchy initial problem with such system it's alright. The problem is fsove involving ode and bvp.
You know Loren's system? It's a first order with system of 3 equations. You know Navier Stoke? It's a first order with square non-linearity term, and people still strugling to understand it.
Your system is 6th order with a power 4 term inside the equation.
Bruno Luong
on 10 Apr 2022
Edited: Bruno Luong
on 10 Apr 2022
This seems to converge on my PC (R2022a, Intel) however not on TMW server.
I iterate on eta span (in a log fashion) and hope the solution of the previous step is a good staring point.
eta0 = 20;
yinit = [3.3517;0;0.5079;0;0.1194;0];
nbiter = 10;
eta0tab = logspace(0,log10(eta0),nbiter);
for k=1:nbiter
eta = linspace(0,eta0tab(k),100);
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
yinit = sol.y(:,1);;
end
plot(eta,F);
function dF = fun_ode(eta,F)
dF=[F(2:6);
(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
function res = bcfun(ya,yb)
res =[ya([2 4 6]); yb(1:3)];
end
Wissem
on 10 Apr 2022
Bruno Luong
on 10 Apr 2022
"However, your last plotting provides the good shape (good values ?). "
I'm not sure as you what mean "good shape" and "good value". I don't know the analytic solution, and I don't know how is other 5 derivative values are correctly matched. Again Imposing condition on fifth derivative is just insane. The derivative is very instable operator, take it five/six time and doing calculation on it is hopeless IMO.
"There is a so strong dependance to the choice of "yinit"...
Because of non linearity, I already tell you that.
When more few algorithms tells that the Jacobian is singular, it just indicates that the problem is so illposed and the dependance of the final state (or at least some subspace of it) to the initial state is numerically zero. It migh be due that your problem is badly scale (for large eta0 >> 1 as I state above) or in the intrincic nature of the problem you are trying to solve.
Also, could you please tell me why did you write the residuals this way ?
res =[ya([2 4 6]); yb(1:3)];
From the bc you are telling us, if you miss my comment in the first code here again
res =[ya(2); % F'(0)=0
ya(4); % F'''(0)=0
ya(6); % F'''''(0)=0
yb(1:3)]; % F(eta0)=0, F'(eta0)=0, F''(eta0)=0,
Bruno Luong
on 10 Apr 2022
I'm curious to know how this equation of sixth order derivative is derived.
To my encounter in physics I rarele see anything more than second order.
@Bruno Luong This is deduced from an equation describing the height field in a so-called "lubrification theory" : the high order is because we consider a pressure component that is of order 4.
The physical problem involves 6 BC on the interval 
that are :
that are : 
@Torsten, The previous code gives a shape that is similar than the one expected but I don't see how I can add new constraint to my problem in order to reduce this dependancy between the final result and the initial guess or the value of η.
% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
eta0 = 40;
yinit = [1;0;0.5079;0;0.1194;0]; % Initial conditions for the bvpinit solver.
nbiter = 10; % Number of iterations.
% We choose to iterate in a log scale the values of eta
% We integrate by using the bvp4c solver for the corresponding ODE's
eta0tab = logspace(0,log10(eta0),nbiter);
for k=1:nbiter
eta = linspace(0,eta0tab(k),100);
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
yinit = sol.y(:,1);
end
% We plot the result
plot(eta,-F);
%We construct a vector of 1-st order ODE's
function dF = fun_ode(eta,F)
dF=[F(2:6);
(1/9*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3];
end
% We built residuals of the method with y(b)-a = 0.
function res = bcfun(ya,yb)
res =[ya([2 4 6]); yb(1:3)];
end
Thank you for your help,
Torsten
on 11 Apr 2022
This is deduced from an equation describing the height field in a so-called "lubrification theory" : the high order is because we consider a pressure component that is of order 4.
Could you give a bibliographic link ?
The previous code gives a shape that is similar than the one expected but I don't see how I can add new constraint to my problem in order to reduce this dependancy between the final result and the initial guess or the value of η.
So you get different converged "solutions" for fixed eta0, but different initial guesses ?
Wissem
on 12 Apr 2022
Bruno Luong
on 12 Apr 2022
I don't understand, the equation that looks similar what your are trying to solve is in the paragraph under equation (4). It is already normalized by the scale (so my comment on the scale is already applied), but I did not recognize exactly the equation or the bc.
Wissem
on 12 Apr 2022
@Bruno Luong The author of this article suggested a form of the solution that is not similar to mine. I have operated the variable changement : 
with 
. Thus, by injecting this form to the previous complete (Eq.1) in a cartesian framework we obtain the original equation of my post. Anyway, BC of F are those verified by h.
with . Thus, by injecting this form to the previous complete (Eq.1) in a cartesian framework we obtain the original equation of my post. Anyway, BC of F are those verified by h.
Bruno Luong
on 12 Apr 2022
So you perform some kind of mapping from infinity domain to finite one with some power law?
By doing so you might render the original scalable problem into a singular Jacobian, and make the numerical solution impossible to solve, IMHO.
Wissem
on 12 Apr 2022
Bruno Luong
on 12 Apr 2022
It appears in the ode as you state in your question
Warning: Matrix is singular to working precision.
> In newtonRaphson2 (line 23)
In shoot4nl (line 14)
as well as under MATLAB
sol = bvp5c(@(xi,G) fun_ode(xi,G,etascale), @bcfun, solinit);
Error using bvp5c
Unable to solve the collocation equations -- a singular Jacobian encountered
Wissem
on 12 Apr 2022
Bruno Luong
on 12 Apr 2022
Why you are certain? Do you know what the effect of t^alpha and t^-beta? Anyone has studied careful the ode and made sure it is solvable, has unique solution and stable?
MATLAB does lies when it throw out an error/warning message.
Torsten
on 12 Apr 2022
From the physical standpoint, is it senseful to assume that there exists a steady-state for equation (1) independent of the initial conditions ?
But you try to calculate a steady-state of the process. And you write that it's not possible ?
@Torsten It's not exactly a "steady-state" since the new variable imposed contains the time dependancy (
). I think there is a little confusion with the words used : Personnally, I prefer to call it a "similarity" solution for this particular formulation.
). I think there is a little confusion with the words used : Personnally, I prefer to call it a "similarity" solution for this particular formulation.
Bruno Luong
on 12 Apr 2022
Edited: Bruno Luong
on 12 Apr 2022
@Wissem-Eddine KHATLA But why do you prefer remarametrization of big-bang type t^-beta*x and big-chill scale t^alpha instead of traveling wave xi = (x - c*t)/Lp as the author of the paper suggests ?
Bruno Luong
on 13 Apr 2022
Edited: Bruno Luong
on 13 Apr 2022
@Wissem-Eddine KHATLA "I don't see how I can add new constraint to my problem in order to reduce this dependancy between the final result and the initial guess or the value of η."
then there is a snip of code that I'm not quite sure what it supposes to demonstrate. On my side I obtain this result
eta0 = 40;
Finit = [1;0;0.5079;0;0.1194;0]
...
Something I don't get is that the author of the paper imposes bc F=1 for xi -> infinity, meaning the steady state where the solution converges is the constant thickness h0, which make sense physically. In your case there is no such thing and solutions can be positive or negative (bc are all 0s), often time numerical solution returns negative solution, and as Torsen has noticed F=0 is a solution as well. I'm not sure how you interpret it with your odd and strong non-linear reparametrization x*t^-beta.
Wissem
on 13 Apr 2022
Bruno Luong
on 13 Apr 2022
Your system seems to have multiple solutions and not stable. IMO it's not a numerical issue but it's just ill-posed problem as it's stated. Sometime such system is derived from energy minimma equilibrium which is lost when one consider the equation alone. May be you should come back to the energey formulation.
Wissem
on 15 Apr 2022
@Bruno Luong , @Torsten : Thank you for your contributions. I came up with a new idea and and a slight modification on the original equation. I have introduced :
- an experimental constant A
- a constraint imposed on the volume that gives a new equation :
Both A and Q are known. Therefore, I obtain these two equations to implement as an ODE's system :
with : 
I have simply added this second equation into the ODE's system in the last line
dF=[F(2:6);
((1/(9*A))*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3; F(1)];
The boundary conditions are still the same :
I just add the condition on the integral that I write : 
I think that allos to have a well posed problem now.
If I test the modified code with : A = Q = 1 : I obtain a plotting. However, if I choose to apply the real values of A and Q : the code crashes again...
%% Solving the problem of elastic droplet's spreading in an elastic regime in the [0;eta0] interval
% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
% Constraint on the volume
eta0 = 50;
yinit = [1; 0; 0.5079; 0; 0.1194; 0; 1]; % Initial conditions for the bvpinit solver.
%yinit = [0.25;0;1;0;0.5;0;1];
nbiter = 10; % Number of iterations.
% We choose to iterate in a log scale the values of eta
% We integrate by using the bvp4c solver for the corresponding ODE's
eta0tab = logspace(0,log10(eta0),nbiter);
for k=1:nbiter
eta = linspace(0,eta0tab(k),100);
solinit = bvpinit(eta,yinit);
sol = bvp5c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
yinit = sol.y(:,1);
end
% We plot the result [F; F'] in a new figure
figure()
plot(eta,-F,'-','LineWidth',2);
grid on
hold on
plot(eta,sol.y(2,:),'LineWidth',2)
legend('F',"F'",'Location','northeast')
xlabel('\eta','FontSize',13);
ylabel("F(\eta),F'(\eta)",'fontsize',13)
% We first construct a vector of 1-st order ODE's
% The last one is the constraint term
function dF = fun_ode(eta,F)
A = 1; %3.7926e-7; % Experimental constant
dF=[F(2:6);
((1/(9*A))*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3; F(1)];
end
% We built residuals of the method with y(b)-a = 0.
function res = bcfun(ya,yb)
Q = 1; %0.4e-9; % Experimental flow rate (m^3/s)
res =[ya([2 4 6]); yb(1:3); yb(7)-Q];
end
Do you guys have an idea about this issue ? Have I correctly written my modifications ? Or is it a problem with my formulation again ?
Bruno Luong
on 15 Apr 2022
You need to figure out the suitable rescaling and work with the equivalent unitless system and not with physics SI unity, which can be result a badly scaled system and implies numerical problem.
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
This code doesn't crash with your experimental data, however there is a warning. In fact your system id 6th order (fake 7th order) and you want to impose 6bc + 1 integral constraint, it cannot meet all in general.
clear
close all
%% Solving the problem of elastic droplet's spreading in an elastic regime in the [0;eta0] interval
% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
% Constraint on the volume
eta0 = 50;
A = 3.7926e-7;
Q = 0.4e-9;
%Finit = [1; 0; 0.5079; 0; 0.1194; 0; 1]; % Initial conditions for the bvpinit solver.
Finit = [-1.0513e+04;0;158.4915;0;-3.6654;0;-1.0513e+04];
nbiter = 20; % Number of iterations.
% We choose to iterate in a log scale the values of eta
% We integrate by using the bvp4c solver for the corresponding ODE's
eta0tab = logspace(log10(1),log10(eta0), nbiter);
for k=1:nbiter
fprintf('iteration %d/%d\n', k, nbiter);
etascale = A^(1/6);
eta = linspace(0,eta0tab(k),100);
xi = eta/etascale;
Ginit = Finit .* etascale.^[0:5 0]';
solinit = bvpinit(xi,Ginit);
try
sol = bvp5c(@(xi, G) fun_ode(xi,G,A,Q,etascale), ...
@(ya,yb) bcfun(ya,yb,A,Q,etascale), ...
solinit);
catch
% if k==nbiter
% fprintf('Fail at the last iteration\n');
% return
% end
continue
end
Finit = sol.y(:,1) ./ etascale.^[0:5 0]';
end
xi = sol.x;
eta = etascale*xi;
F = sol.y(1,:);
Fp = sol.y(2,:) / etascale;
% We plot the result [F; F'] in a new figure
figure()
plot(eta,-F,'-','LineWidth',2);
grid on
hold on
plot(eta,Fp,'LineWidth',2)
legend('F',"F'",'Location','northeast')
xlabel('\eta','FontSize',13);
ylabel("F(\eta),F'(\eta)",'fontsize',13)
% We first construct a vector of 1-st order ODE's
% The last one is the constraint term
function dG = fun_ode(xi,G,A,Q,etascale)
eta = xi*etascale;
F = G ./ (etascale).^[0:5 0]';
dF=[F(2:6);
((1/(9*A))*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3;
F(1)/A];
dG = dF .* (etascale).^[1:6 1]';
end
% We built residuals of the method with y(b)-a = 0.
function res = bcfun(ya,yb,A,Q,etascale)
ya = ya ./ (etascale).^[0:5 0]';
yb = yb ./ (etascale).^[0:5 0]';
res =[ya([2 4 6]); yb(1:3); yb(7)-Q/A];
end
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
@Wissem-Eddine KHATLA also it seems I(0) must be 0, so that I(eta0) = Q, but I don't see where you impose I(0)=0 in your code. May be what you should impose is
I(eta)-I(0) = Q
Torsten
on 15 Apr 2022
I don't know if bvp4c accepts this destruction of the banded structure of the internal Jacobian, but yes, the last condition should be
yb(7) - Q/A - ya(7)
And starting should be better with
Finit = [-1.0513e+04;0;158.4915;0;-3.6654;0;0];
instead of
Finit = [-1.0513e+04;0;158.4915;0;-3.6654;0;-1.0513e+04];
I guess.
Bruno Luong
on 15 Apr 2022
There is something very wrong with the so called experimental data A and Q provided.
@Wissem-Eddine KHATLA need to go back to the blackboard and check the model/parameter carefully.
I'm also surprised H integral doesn't have eta^(n-1) term (volume/mass conservation) as if it's 1D problem.
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
@Torsten: "I don't know if bvp4c accepts this destruction of the banded structure of the internal Jacobian, but yes, the last condition should be"
yb(7) - Q/A - ya(7)
No it doesn't work, bvpxc detects { yb(7)-ya(7) = 0 } is a null space of the Jacobian and throws an error. In short one cannot add integral constraint by adding a dummy variable.
Wissem
on 15 Apr 2022
@Bruno Luong Thank you for your remarks. In fact, by following your advice, we can write a unitless equation in which A = Q = 1. And thus perform my old code that gives this :
also it seems I(0) must be 0, so that I(eta0) = Q, but I don't see where you impose I(0)=0 in your code. May be what you should impose is
I have thaught about it : but I have a 6th-order equation + an integral "constraint" so I thaught that 7 BC should be enough for the solver to run isn't it ?
need to go back to the blackboard and check the model/parameter carefully.
I'm also surprised H integral doesn't have eta^(n-1) term (volume/mass conservation) as if it's 1D problem.
I am not sure to understand why the integral should have eta(n-1) terms as the constraint should be applied on the whole interval.
Thank you in advance,
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
I don't know what mean eta, it looks like a radial spatial variable, and when one integrate in cylindrical/spherical coordinate one need to put r ^(n-1) inside the integral, where n is the dimension of the space.
If you don't impose I(0)=0, you don't actually impose integral constraint, since bvpxc solvers simply changes freely I(0) so as to match I(eta0) to Q. Mathematically you don't impose anything new in your equation. Somehow bvpxc behaves better, but IMO it's just a pure luck.
Wissem
on 15 Apr 2022
Bruno Luong
on 15 Apr 2022
Try it, it won't work, see my previous answer to Torsen.
Wissem
on 15 Apr 2022
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
I believe you cannot, which is unfortunate but actually quite logic. You cannot add a constraint to a solution of the ode with all the bc prescribed, even if the solution looks unstable (and if I believe in the math paper you posted earlier, the system is unstable).
I remember we started with the domain of integration [-eta0,eta0] and the conditions
F(-eta0) = F'(-eta0) = F(eta0) = F'(eta0) = F'(0) = F'''(0) = 0.
Then we changed the domain to [0,eta0] assuming symmetry of the solution. In my opinion, this only uses
F'(0) = F'''(0) = F'''''(0) = F(eta0) = F'(eta0) = 0.
Since I don't know where the extra condition
F''(eta0) = 0
stems from, it might be possible to drop this condition, add the equation dF(7)/d(eta) = F(1) to the system and add the two boundary conditions
ya(7) = 0, yb(7) - Q = 0
Just a suggestion.
Wissem
on 15 Apr 2022
Bruno Luong
on 15 Apr 2022
Edited: Bruno Luong
on 15 Apr 2022
If you want to add an (integral) strict equality, you either have to remove one bc, or add a free parameter. You simply cannot add a condition to a set of discrete solutions fully defined by the original ode. This is just mathematically wrong.
Wissem
on 15 Apr 2022
Bruno Luong
on 15 Apr 2022
"If I have correctly understood : the last value in the vector
yinit = [1; 0; 0.5079; 0; 0.1194; 0; 1]; % Initial conditions for the bvpinit solver.
is the guess for the integral value isn't it ?"
Not exactly, since yinit will be replicate on eta span, so it is a guess of both 0 and the integral.
Wissem
on 15 Apr 2022
Bruno Luong
on 15 Apr 2022
I simply state the doc of bvpinit
"Vector – For each component of the solution, bvpinit replicates the corresponding element of the vector as a constant guess across all mesh points. That is, yinit(i) is a constant guess for the ith component yinit(i,:) of the solution at all the mesh points in x."
Since the boundary conditions will be
res =[ya([2 4 6 7]); yb(1:2); yb(7)-Q];
it would be reasonable to set the initial guess for y(7) to
Q*eta/eta0
because this ensures that the initial guess is consistent with the boundary conditions imposed.
Instead of giving a constant vector yinit, the conditions now had to be specified in a function:
solinit = bvpinit(eta,@(eta)guess(eta,eta0,Q);
function g = guess(eta,eta0,Q)
g = [1; 0;0.5079; 0; 0.1194; 0; Q*eta/eta0];
end
Wissem
on 16 Apr 2022
Bruno Luong
on 16 Apr 2022
Edited: Bruno Luong
on 16 Apr 2022
@Wissem-Eddine KHATLA This is not correct anonymous function syntax.
@guess(eta,eta0,Q)
It looks like you want to write (without @)
bvpinit(eta, guess(eta,eta0,Q))
which is vector initialization syntax.
IMO Torsen wants to use function for bcpinit since i will not replicate value on the eta-span: from the doc
- Function – For a given mesh point, the guess function must return a vector whose elements are guesses for the corresponding components of the solution. The function must be of the form y = guess(x)x is a mesh point and y is a vector whose length is the same as the number of components in the solution. For example, if yinit is a function, then at each mesh point bvpinit callsy(:,j) = guess(x(j))For multipoint boundary value problems, the guess function must be of the formy = guess(x, k)y is an initial guess for the solution at x in region k. The function must accept the input argument k, which is provided for flexibility in writing the guess function. However, the function is not required to use k.
I have tried all Torsen's ideas not not post here since none of them bring anything, not saying it degrade the stability.
After playing quite a bit of it IMO forcing the integral is a false good idea.
Wissem
on 16 Apr 2022
@Bruno Luong Thank you for this post. To be honest, I don't see any other alternative to the forcing integral to solve this problem. I am thinking about adding an unknown parameter which would physically corresponds to the initial height 
and with the parametrization : 
. This would imply 2 changements in the set of equations :
and with the parametrization : . This would imply 2 changements in the set of equations :
This changement implies that we would need a new bc to "close" the problem. Maybe like this ?
@Bruno Luong Do you think it would made a big difference to write it this way ?
Thank you again for your useful advices,
Bruno Luong
on 16 Apr 2022
Edited: Bruno Luong
on 16 Apr 2022
I think the integral is overly constrained the problem mathematically, then you mix experimental data with theoretical model with 0 error, this will not work since you always have errors in both sides (your Q value of 1e-9 is totally off IMO).
The authors of the paper seem to study throughly the equation and arrive to compute solution with MATLAB, why don't you follow exactly their approach?
Wissem
on 16 Apr 2022
@Bruno Luong It is because we have not the same experimental approach. As you may have noticed, they build their solution in order to "match" it to a film of thickness 
(this notation is different from the one of my previous post in which I've noted 
) : something that I don't have here (
). I am not sure about how their parametrization should work in this case
(this notation is different from the one of my previous post in which I've noted ) : something that I don't have here (). I am not sure about how their parametrization should work in this case
Torsten
on 16 Apr 2022
solinit = bvpinit(eta,@(eta)guess(eta,eta0,Q));
I mean : what's the difference between this line and this one ??
solinit = bvpinit(eta,@guess(eta,eta0,Q));
The difference is that the first line will work with bvp4c, the other will not.
Did you try ?
If you define the initial guess for the solution of your ODE by a function, it has the form
function Finit = guess(eta)
...
end
Now you changed the input list to
function Finit = guess(eta,eta0,Q)
In order to tell MATLAB at which position of the list of inputs it should pass eta, it is necessary to indicate this as
@(eta) guess(eta,eta0,Q)
Torsten
on 16 Apr 2022
It is because we have not the same experimental approach. As you may have noticed, they build their solution in order to "match" it to a film of thickness (this notation is different from the one of my previous post in which I've noted ) : something that I don't have here (). I am not sure about how their parametrization should work in this case
Although the underlying physical problem in the article might be different, I'd just try whether the MATLAB approach with bvp4c works for their problem. This would at least be a confirmation that the problems you encounter at the moment are due to your problem, not due to MATLAB.
@Torsten I was about to do it but I notice that they are solving an equation of order 5 with 4 boundary conditions (the one noted 
). Moreover, they seem to have a "hidden" unknown parameter nammed A and they write that its value is 1.35 after numerical calculation... So they must have 6 BC but only 4 are detailed
). Moreover, they seem to have a "hidden" unknown parameter nammed A and they write that its value is 1.35 after numerical calculation... So they must have 6 BC but only 4 are detailed
Torsten
on 16 Apr 2022
Then I'd contact the authors to get more infos.
Bruno Luong
on 16 Apr 2022
They claim using bvp4c.
@Torsten @Bruno Luong I've came up with news about the paramtrization shown in the article. I am aware that is certainly a different problem than the one exposed in my original post but since we are discussing it from the beginning here is where I have some issues :
The differential problem solved by the authors is : 
.
.3 boundary conditions are easily exposed :
I am wondering how I can implement the 2 other ones since they are deduced after a new parametrization by writing : 
with 
. Asymptotically, ϕ verifies the equation : 
which have 5 exponential solutions of the form 
and 
. It implies that 2 other conditions are deduced for F using ϕ. ϕ should verify 2 differential equations that the authors consider as boundary conditions on 
which are :
with . Asymptotically, ϕ verifies the equation : which have 5 exponential solutions of the form and . It implies that 2 other conditions are deduced for F using ϕ. ϕ should verify 2 differential equations that the authors consider as boundary conditions on which are :
I am not sure how a boundary condition solution to an other differential equation could be implemented on MATLAB. The process should be the same :
- Creating a system of 1st order ODE on F.
- Implementing the boundary conditions and residuals on the 3 obvious given BC.
- But how to implement the 2 other ones as a solution to an other differential equation ??
Thank you in advance for your help.
Torsten
on 18 Apr 2022
What is "alpha" ?
Bruno Luong
on 18 Apr 2022
@Torsten alpha is blowup parameter in Evan paper.
@Wissem-Eddine KHATLA Sorry, I give up, I'm not expert in this physics of peeling, viscous, lubrification model to be able to help you.
Wissem
on 18 Apr 2022
I can only naively evaluate the two expressions you give ( I did not look into the article):
(D+1)(D-alpha)(D+alpha') phi =
(D^3+D^2-abs(alpha)^2*D-abs(alpha)^2) phi =
phi''' + phi'' - abs(alpha)^2*phi' - abs(alpha)^2*phi = 0
D(D+1)(D-alpha)(D+alpha') =
phi'''' + phi''' - abs(alpha)^2*phi'' - abs(alpha)^2*phi' = 0
Since F = 1 + phi,
F''' + F'' - abs(alpha)^2*F' - abs(alpha)^2*(F-1) = 0
F'''' + F''' - abs(alpha)^2*F'' - abs(alpha)^2*F' = 0
Since F''' = F'''' = 0 at oo,
F'' - abs(alpha)^2*F' - abs(alpha)^2*(F-1) = 0
- abs(alpha)^2*F'' - abs(alpha)^2*F' = 0
thus
F' + F'' = 0
F'*(1+abs(alpha)^2) + abs(alpha)^2*F - abs(alpha)^2 = 0
Wissem
on 18 Apr 2022
Torsten
on 18 Apr 2022
Why mixed ?
It's just
yb(2)+yb(3)
yb(2)*(1+abs(alpha)^2)+abs(alpha)^2*yb(1)-abs(alpha)^2
What's the problem ?
Wissem
on 18 Apr 2022
@Torsten : It seems like the boundary conditions are the following :
and 
for 
I have tried to implement these conditions for the BVP solver like this :
function res = bc(FL,FR)
alpha = cos(4*pi/5) + i*sin(4*pi/5);
beta = conj(alpha);
gamma = abs(alpha);
res = [FL(4,1) ; FL(5,1);...
FR(1,1) - 1 ;...
FR(1,1)-FL(1,2) ; FR(2,1)-FL(2,2) ; FR(3,1)-FL(3,2) ; FR(4,1)-FL(4,2) ; FR(5,1)-FL(5,2) ;...
(1-(beta+alpha))*FR(3,2) + (-(alpha+beta)+gamma^2)*FR(2,2) ; (-(beta+alpha)+gamma^2)*FR(3,2) + (FR(2,2)-1)*gamma^2];
end
I am just asking for a review on my "res" function : I am a bit disturbed by the the condition : F(0) = 1. Do you think I am writing it correctly ? I saw in documentation that there is a "region" function but I am not sure how to implement in this case.
Thank you in advance @Torsten
Looks ok. Although I don't understand how the F''' and F'''' terms in the boundary condition at eta0 could drop out - now that you don't impose these two conditions at eta0, but at -eta0.
The region parameter is an input parameter to the function fun_ode where you define the differential equation.
With this parameter, you can define different differential equations depending on the region you're in.
In your case, the differential equation for eta in [-eta0,0] could differ from the differential equation in [0,eta0].
You don't need to do anything with this parameter. Nevertheless, it is part of the list of input parameters for fun_ode (eta,F,region).
Wissem
on 19 Apr 2022
Torsten
on 19 Apr 2022
I thought you integrate over [0,oo) and assume symmetry to eta=0 with one boundary condition at eta=0 and four boundary conditions at eta=eta0. This is not true for the problem treated in the paper ?
@Torsten Yes : it was my case but it seems like they integrate on 
in the paper with the conditions that I have detailed previously. Unfortunately, my code doesn't work since I am supposed to find 
...
in the paper with the conditions that I have detailed previously. Unfortunately, my code doesn't work since I am supposed to find ...The code is the following :
eta0 = 20;
etamesh = [linspace(-eta0,0,100),linspace(0,eta0,100)];
Finit = [1; 0.5; 1; 0; 0];
solinit = bvpinit(etamesh,Finit)
sol = bvp5c(@fun_ode, @bc, solinit);
figure()
plot(sol.x,sol.y(1,:),'linewidth',1.5) % Plotting F
ylim([-2 2]);
hold on
plot(sol.x,sol.y(3,:),'linewidth',1.5) % Plotting F''
plot(sol.x,sol.y(5,:),'linewidth',1.5) % Plotting F''''
grid on
legend("F","F''","F''''",'location','northwest')
function dFdeta = fun_ode(eta,F,region)
dFdeta=[F(2);F(3);F(4);F(5);(1-F(1))/F(1)^3];
end
% Boundary conditions :
% F'''(-eta0) = F''''(-eta0) = 0
% F(0) = 1
% [1] = [2] = 0 on eta = eta0
function res = bc(FL,FR)
alpha = cos(4*pi/5) + i*sin(4*pi/5);
beta = conj(alpha);
gamma = abs(alpha);
res = [FL(4,1) ; FL(5,1);...
FR(1,1) - 1 ;...
FR(1,1)-FL(1,2) ; FR(2,1)-FL(2,2) ; FR(3,1)-FL(3,2) ; FR(4,1)-FL(4,2) ; FR(5,1)-FL(5,2) ;...
(1-(beta+alpha))*FR(3,2) + (-(alpha+beta)+gamma^2)*FR(2,2) ; (-(beta+alpha)+gamma^2)*FR(3,2) + (FR(2,2)-1)*gamma^2];
end
I think that everything is alright with my code : maybe you will find an error on it ?
Thank you.
Torsten
on 20 Apr 2022
Looks ok for me (in terms of content, I can't contribute here).
What do you mean by
Unfortunately, my code doesn't work since I am supposed to find F''(-eta0) = 1.35...
Does your code not work or can't you extract F''(-eta0) or does F''(-eta0) not equal 1.35... ?
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
on 20 Apr 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
