**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

54 views (last 30 days)

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 :

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.

##### 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-Eddine KHATLA
on 7 Apr 2022

### Answers (2)

Torsten
on 6 Apr 2022

Edited: Torsten
on 6 Apr 2022

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-Eddine KHATLA
on 6 Apr 2022

@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
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-Eddine KHATLA
on 6 Apr 2022

@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
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-Eddine KHATLA
on 6 Apr 2022

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

Torsten
on 6 Apr 2022

Edited: Torsten
on 6 Apr 2022

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-Eddine KHATLA
on 7 Apr 2022

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-Eddine KHATLA
on 7 Apr 2022

Edited: Wissem-Eddine KHATLA
on 7 Apr 2022

Torsten
on 7 Apr 2022

Edited: Torsten
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-Eddine KHATLA
on 7 Apr 2022

@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
on 7 Apr 2022

Torsten
on 7 Apr 2022

Edited: Torsten
on 7 Apr 2022

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-Eddine KHATLA
on 7 Apr 2022

@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
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-Eddine KHATLA
on 8 Apr 2022

@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
on 8 Apr 2022

Edited: Torsten
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-Eddine KHATLA
on 8 Apr 2022

@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-Eddine KHATLA
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-Eddine KHATLA
on 8 Apr 2022

Torsten
on 8 Apr 2022

Edited: Torsten
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 ?

Torsten
on 9 Apr 2022

Edited: Torsten
on 9 Apr 2022

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-Eddine KHATLA
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 :

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
on 9 Apr 2022

Edited: Torsten
on 9 Apr 2022

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-Eddine KHATLA
on 9 Apr 2022

Edited: Wissem-Eddine KHATLA
on 9 Apr 2022

@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
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-Eddine KHATLA
on 9 Apr 2022

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

Wissem-Eddine KHATLA
on 10 Apr 2022

Edited: Wissem-Eddine KHATLA
on 10 Apr 2022

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

Bruno Luong
on 9 Apr 2022

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

Wissem-Eddine KHATLA
on 10 Apr 2022

Edited: Wissem-Eddine KHATLA
on 10 Apr 2022

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

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

Wissem-Eddine KHATLA
on 10 Apr 2022

Edited: Wissem-Eddine KHATLA
on 10 Apr 2022

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-Eddine KHATLA
on 10 Apr 2022

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

Wissem-Eddine KHATLA
on 11 Apr 2022

Edited: Wissem-Eddine KHATLA
on 11 Apr 2022

@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
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-Eddine KHATLA
on 12 Apr 2022

Bruno Luong
on 12 Apr 2022

Wissem-Eddine KHATLA
on 12 Apr 2022

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-Eddine KHATLA
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-Eddine KHATLA
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

Wissem-Eddine KHATLA
on 12 Apr 2022

Edited: Wissem-Eddine KHATLA
on 12 Apr 2022

Wissem-Eddine KHATLA
on 12 Apr 2022

Edited: Wissem-Eddine KHATLA
on 12 Apr 2022

Bruno Luong
on 12 Apr 2022

Edited: Bruno Luong
on 12 Apr 2022

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-Eddine KHATLA
on 13 Apr 2022

@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
on 13 Apr 2022

Wissem-Eddine KHATLA
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

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-Eddine KHATLA
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-Eddine KHATLA
on 15 Apr 2022

@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

Wissem-Eddine KHATLA
on 15 Apr 2022

Bruno Luong
on 15 Apr 2022

Edited: Bruno Luong
on 15 Apr 2022

Torsten
on 15 Apr 2022

Edited: Torsten
on 15 Apr 2022

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-Eddine KHATLA
on 15 Apr 2022

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

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

Thank you,

Bruno Luong
on 15 Apr 2022

Edited: Bruno Luong
on 15 Apr 2022

Wissem-Eddine KHATLA
on 15 Apr 2022

@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
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-Eddine KHATLA
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."

Torsten
on 15 Apr 2022

Edited: Torsten
on 15 Apr 2022

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-Eddine KHATLA
on 16 Apr 2022

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

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-Eddine KHATLA
on 16 Apr 2022

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.

Wissem-Eddine KHATLA
on 16 Apr 2022

Edited: Wissem-Eddine KHATLA
on 16 Apr 2022

Wissem-Eddine KHATLA
on 18 Apr 2022

Edited: Wissem-Eddine KHATLA
on 18 Apr 2022

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

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-Eddine KHATLA
on 18 Apr 2022

Edited: Wissem-Eddine KHATLA
on 18 Apr 2022

Wissem-Eddine KHATLA
on 18 Apr 2022

Torsten
on 18 Apr 2022

Edited: Torsten
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-Eddine KHATLA
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-Eddine KHATLA
on 18 Apr 2022

Wissem-Eddine KHATLA
on 19 Apr 2022

Edited: Wissem-Eddine KHATLA
on 19 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

Torsten
on 19 Apr 2022

Edited: Torsten
on 19 Apr 2022

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-Eddine KHATLA
on 19 Apr 2022

Torsten
on 19 Apr 2022

Wissem-Eddine KHATLA
on 20 Apr 2022

Edited: Wissem-Eddine KHATLA
on 20 Apr 2022

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

Wissem-Eddine KHATLA
on 20 Apr 2022

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

Compared to the one of the article :

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)