Solving a 6th order non-linear differential equation in Matlab

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 :

with :

According to the boundary conditions we thus have : 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...

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.

[Bruno Luong]

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]

@Bruno Luong Thanks for your reply. @Torsten provided a correction to my first attempt by using various in-built functions but it seems like I got an error stating that the "fsolve" function is not correct...

Answer 1

[Torsten]

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

[Wissem]

@Torsten Thank you again for your help. I've tried to use the bvp4c in-built function but didn't succeed so I came back to my own functions...

Unfortunately, I didn't understand the aim of your function "res" ? Because "etaspan" is not built as I expected to be for the ode45 solver ? Also : what is the objective of defining "bc" like this ?

Thank you in advance.

Best regards.

[Torsten]

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]

@Torsten How do you think we can adjust the step size of integration in this case ? Because, while running your script, I obtain the following error :

Warning: Failure at t=-9.999884e+00.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (2.842171e-14) at time t. 

For instance, I used to do it manually in my previous script.

Thanks

[Torsten]

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]

I wanted to adjust the stepsize precisely because I got an error on it on Matlab :

Warning: Failure at t=-9.999884e+00.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (2.842171e-14) at time t. 

It seems like the solver is not able to adjust it efficiently maybe ? Or maybe is it a difference between the two solvers in Matlab and Octave.

[Torsten]

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.

[Torsten]

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]

Hi @Torsten : Thank you for your great contribution. Unfortunately, for some reason, I am not able to obtain the same plot as yours. I got an error that states :

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.

<stopping criteria details>

I first tried your script directly but then I changed some values on bc0 : it didn't provide any result...

Could you please take the time to explain the modifications that you provided from your previous correction ?

Thank you in advance,

Best regards.

[Wissem]

@Torsten I still have some issues with the "fsolve" in-built function and don't know what's the reason behind... I am still trying to figure it out.

[Torsten]

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]

@Torsten You mean it would be more convenient to change the variable eta into something else in order to avoid the +/-inf condition ?

Also, the shape that you obtain is the one that is expected. I don't understand why I keep facing this error about "fsolve" I tried to run on octave and I have an error too. There is something that I am missing.

[Torsten]

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.

[Torsten]

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]

@Torsten I am currently using the version 6.3.0. More precisely, I am getting this error :

error: invalid call to script /Users/wissem-eddine/Documents/MATLAB/fun.m
error: called from
    fun
    /private/var/folders/4g/s18zdpyn2s5bfgr15s40n21c0000gn/T/octave/octave_VINHsb.m>@<anonymous> at line 3 column 18
    fsolve at line 242 column 8

[Torsten]

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]

@Torsten Thank you for your contribution on the bvp4c script : I was wondering if it was easier to use it for my purpose. However, I would like to fully understand the first script since I keep encountering an error on "fsolve":

• Why did you choose to compute 4 residuals in the function ?

function res = fun(bc,etaspan)

Since we "only" have to check for F=F'=0 at eta0 (and -eta0 after)

• 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 ?

Thank you again for sharing and for your help,

Best regards.

[Torsten]

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]

@Torsten Thank you for your explanations. I tend to forget that my problem is somehow a "3 points problem"...

I tried your "bvp4c" script : it gives an error stating that a singular matrix is encountered :

Error using bvp4c (line 197)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in shoot7nl (line 5)
sol = bvp4c(@fun_ode, @bc, solinit);

It is surely due to the initial conditions values provided I think.

[Wissem]

@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.

[Torsten]

fsolve is not equal fsolve - the solvers are different.

Octave only took the same names as in MATLAB to indicate that both solvers handle the same sort of problems and in order to be compatible with MATLAB code.

The same is true for ODE45.

[Wissem]

@Torsten Do you have some advice to avoid this error ? Because I am still trying to solve it by modifying etaspan and b0 but I haven't succeed yet. I want to fully understand it for a MATLAB use. Thank you again.

[Torsten]

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 ?

[Torsten]

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]

@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 :

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...

Thank you again.

[Torsten]

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

[Wissem]

@Torsten I would also like to ask you how is the residual function "bc" built in your "bvp4c" script please ? More precisely : I don't see why we need to precise two inputs (FL and FR) ?

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   

P.S : Do you have a "bvpinit" package for OCTAVE to advise ?

Thank you again.

[Torsten]

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]):

https://de.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

"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]

@Torsten Thank you for the reference : Indeed, this case was not detailed in the documentary section of MATLAB. Also, have you ever encountered this error with this solver ?

Error using bvp4c (line 197)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Thank you,

[Torsten]

https://de.mathworks.com/matlabcentral/answers/120468-why-do-i-receive-an-error-about-a-singular-jacobian-when-i-use-the-bvp4c-function-within-matlab

I suggest you start with the solution obtained from the shooting method as initial guess.

I'm curious whether bvp4c will reproduce the solution.

[Wissem]

@Torsten I tried the bvp4c solver using the shooting method result : it doesn't seem to work well with the following code :

% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0

eta0 = 20;
eta = linspace(0,eta0,100);

yinit = [3.3517;0;0.5079;0;0.1194;0]; % solution of the OCTAVE ode45 integration
solinit = bvpinit(eta,yinit);
sol = bvp4c(@fun_ode, @bcfun, solinit);
eta = sol.x;
F = sol.y(1,:);
figure
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); ya(4); ya(6); yb(1); yb(2); yb(3)];
end

Another user suggests me that the problem is not stiff : I am confused since a similar equation is treated in the litterature so I think this is still possible for me to solve it.

Thank you,

Answer 2

[Bruno Luong]

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

[Torsten]

@Bruno Luong

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),

[Wissem]

@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.

Thank you again,

[Torsten]

@Wissem-Eddine KHATLA

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]

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]

@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.

@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]

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).

[Wissem]

@Bruno Luong You are right about the "insane" part of this problem. In my original post, I tried a home-made function in order to solve it but thus @Torsten helped me to use in-built functions of MATLAB : it produced a bell-shape plotting which is quiet the shape expected. However, there is a stability problem using "ode45" solvers that's why I was advised to use "bvp4c" : I am surprised to see that "ode45" and "bvp4c" are not giving the same solution at all.

[Bruno Luong]

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]

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]

@Bruno Luong I think you are probably right about the stifness of this problem... However, your last plotting provides the good shape (good values ?). There is a so strong dependance to the choice of "yinit"...

Also, could you please tell me why did you write the residuals this way ?

res =[ya([2 4 6]); yb(1:3)];

Thank you again,

[Bruno Luong]

"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,

[Torsten]

I think one main problem in the solution of this BVP is that F = 0 is also a solution apart from the nontrivial one we are looking for.

Maybe there are infinitely many ... who knows ?

[Bruno Luong]

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.

[Wissem]

@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 :

• 
• 

@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]

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]

@Torsten Sure : you could find the article attached to this post. The equation from which the problem is derived is the (Eq.1) in the second page with the boundary conditions details. I hope that it could help to understand the mathematical formalism of my problem.

[Bruno Luong]

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]

@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.

[Bruno Luong]

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]

@Bruno Luong I am not sure to see how the original problem is "different" from the one that I wrote : where does the singular Jacobian appear ?

[Bruno Luong]

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]

@Bruno Luong Yes but I mean there is no "additional" reason that could justify the apparition of this singular matrix in my problem rather than the original one

[Bruno Luong]

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]

@Wissem-Eddine KHATLA

From the physical standpoint, is it senseful to assume that there exists a steady-state for equation (1) independent of the initial conditions ?

[Wissem]

@Torsten Unfortunately, in this case : it's not possible... But it would have made the problem easier to write indeed.

[Torsten]

@Wissem-Eddine KHATLA

But you try to calculate a steady-state of the process. And you write that it's not possible ?

[Wissem]

@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.

[Bruno Luong]

@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]

@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]

@Bruno Luong The fact is that they impose F=1 in infinity because the solution converges to a pre-existed thickness h0 : something that is not true in my configuration. I have conditions of zero-contact angle and zero curvature at the interface η which makes it a bit different...

I am just wondering if this configuration is possible numerically : I don't know tbh. I have found an article that treats mathematically about this kind of problems but I haven't read it yet.

[Bruno Luong]

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]

@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 ?

Thank you again for your help @Torsten and @Bruno Luong

[Bruno Luong]

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]

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]

@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]

@Wissem-Eddine KHATLA

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]

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]

@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]

@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]

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]

@Bruno Luong Do you think I should write the condition I(0) = 0 like this :

yb(7) - ya(7) - Q

Is it correct ?

Thank you

[Bruno Luong]

Try it, it won't work, see my previous answer to Torsen.

[Wissem]

@Bruno Luong Sorry, I didn't see it. How do you think I should write it ? Because I've tried unsuccessfully..

[Bruno Luong]

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).

[Torsten]

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]

@Bruno Luong I've found a post on another forum where they did the same

https://scicomp.stackexchange.com/questions/26665/how-to-implement-an-integral-condition-when-solving-a-bvp-in-matlab

I really think that imposing this constraint is important for the physical meaning of the resolution...

Thank you,

[Bruno Luong]

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]

@Bruno Luong , @Torsten I've applied Torsten suggestion and deleted a bc (F''(eta0) = 0). @Torsten I've added this conditions since I was looking for a bc that would "close" the system : since it means that I have no curvature at infinity, I tried F''(eta0) = 0.

I have a doubt on the way that I am supposed to built the vector for initial conditions since the new system of ODE's is :

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

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 ?

Thank you,

[Bruno Luong]

"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]

@Bruno Luong I am not sure to understand how is it possible to be a simultaneous guess for 0 and the integral at the same time ? You mean I(0) and I(eta0) ?

[Bruno Luong]

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."

[Torsten]

@Wissem-Eddine KHATLA

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]

@Torsten I am sorry to insist on this point again but I really don't get the aim of the line :

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));

I understand it the same way but I am unable to perceive the difference even after consulting the related documentation: could you tell me again why do you write it this way ?

Thank you again,

[Bruno Luong]

@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]

@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 :

• 
• 

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]

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]

@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

[Torsten]

@Wissem-Eddine KHATLA

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]

@Wissem-Eddine KHATLA

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.

[Wissem]

@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

[Torsten]

Then I'd contact the authors to get more infos.

[Bruno Luong]

They claim using bvp4c.

[Torsten]

But the boundary conditions used seem unclear. Or did you succeed to extract the mathematical formulation of the problem from the article ?

[Wissem]

@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 :

• 
• 

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]

What is "alpha" ?

[Bruno Luong]

@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]

@Torsten Sorry : α is the complex solution of that has a negative real part.

There are 5 solutions : the obvious one is -1. 2 solutions have positive real parts and 2 others have negative real parts (α and )

[Wissem]

@Bruno Luong Understandable but my question is a MATLAB request (and my previous post was more about mathematics rather than physics). Anyway, thanks for your help !

[Torsten]

@Wissem-Eddine KHATLA

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]

@Torsten Thank you for this developpement : but my question was about how to implement this kind of condition on MATLAB for a solver like bvp4c ? I mean is it even possible since it's a condition that implies "mixed" derivatives ? The 3 other remaining solutions could be easily coded.

[Torsten]

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]

@Torsten Ok : I thaught it was not possible to write yb(1), yb(2)... in the same residual relation. Thanks : I will try to implement it and compare with our previous attempts.

[Wissem]

@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

[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]

@Torsten I would have liked to impose it on eta0 too but I don't know if it wouldn't make the system over constrained ? (7 boundary conditions for a 5th order equation)..

[Torsten]

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 ?

[Wissem]

@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 ...

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]

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... ?

[Wissem]

@Torsten F''(-eta0) doesnt equal 1.35. For instance, here is my plot :

Compared to the one of the article :