6th order differential equation [bvp5c issue]

Hello everyone,

I am trying to solve a 6-th order non-linear differential equation :

with Fa height function defined on

These are the 6 boundary conditions needed to solve it :

This is similar to my previous post and after several attempts and a very appreciated help, I am using the BVP5C solver in order to integrate the following ODE system such as :

with :

Here is the script of the code :

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

eta0 = 1;
yinit = [1; 0; 0.5079; 0; 0.1194; 0]; % Initial conditions for the bvpinit solver.
nbiter = 45; % 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

plotting = figure(4)
plot(eta,-F,'-','LineWidth',2);
grid on
hold on
plot(eta,sol.y(2,:),'LineWidth',2)
legend('F',"F'",'fontsize',13,'fontname','Times','Location','northeast')
xlabel('\eta','FontSize',13,'FontName','Times');
ylabel("F(\eta), F'(\eta)",'fontsize',13,'fontname','Times')


% We first construct a vector of 1-st order ODE's
% The last one is the constraint term

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

Another users of the forum advised me to use this logarithmic loop in order to ease the integration process.

The problem that I encounter is that : despite imposing a condition , the resulted plot never reach 0 and I keep having some errors about high residuals of the method.

Warning: Unable to meet the tolerance without using more than 1666 mesh points.

Do you think that the solver is not able to find a converged solution despite these boundary solutions ? I have tried to specify manually the general form of the jacobian matrix for bvp5c but it has never worked. I am thinking of using another software but I am not sure if it would be worthy...

Thank you in advance for your help,

Best regards.

Answer 1

[Torsten]

You divide by F(1)^3.

For the function F to be well-behaved at eta = 1, you must have that

5*F(1)-4*F(2)-27*F(1)^2*F(2)*F(6)

has a zero of order at least 3 at eta = 1.

I don't know whether this is solvable.

[Bruno Luong]

Of course it solvable with F==0.

[Torsten]

So here we are again.

[Bruno Luong]

Exactly, before question about failure of numerical method, may be someone must ensure there exists a non-trivial solution.

[Wissem-Eddine KHATLA]

@Bruno Luong I guess if I found an analytical solution to this equation, I wouldn't bother looking for a numerical analysis :) Thanks anyway for your kind advices :)

@Torsten : I am dealing with this singularity from the beginning : I have alreay used the option SingularTerm by the past but never for this kind of issue. As you may remember in my original post, I have provided some "home-made" codes in order to perform the integration manually but it didn't work either. I was hoping that this solver was about to give "at least" a robust solution. Maybe give a try with bvode of Scilab

[Torsten]

Bruno's point was: As long as we don't know whether a solution apart from the trivial solution even exists, there will always be a doubt about "solutions" obtained from a numerical solver, especially because the convergence behaviour was always unstable and because of the singularity at eta = 1.

I tested "bvode" from scilab on one of your earlier problems because "colsys" (its name in fortran times) in my opinion is the best bvp integrator existing. But it didn't come up with a solution, either.

[Wissem-Eddine KHATLA]

@Torsten Yes : thank you for your advice. At least, I think it’s a clearer comment. I am disappointed to hear that you have already tried « bvode » without any success because I was about to test it for this problem.

[Torsten]

Nevertheless, you should give it a try. My attempt was not very thorough.

Answer 2

[Wissem-Eddine KHATLA]

@Torsten : Thank you. I will definitely try the solver : I have nothing to lose by doing so. Also, I do think that there is no analytical solution for this particular equation. I guess that I will be obliged to come back to my original idea of solving it with a finite difference method and resolving the system with a Newton method.

[Torsten]

Also, I do think that there is no analytical solution for this particular equation.

It's not necessary that one is able to explicitly write down an analytical solution. In mathematics, you can often proove existence and uniqueness of solutions for differential equations which have certain properties (e.g. Lipschitz continuity of the function to be integrated) for which nobody is able to write down a solution formula. The authors of the publications you included, e.g., hopefully know that they try to numerically find solutions for problems for which solutions exist.

I guess that I will be obliged to come back to my original idea of solving it with a finite difference method and resolving the system with a Newton method.

But that's what bvp4c does...

[Bruno Luong]

@Torsten: "But that's what bvp4c does..."

What I'm guessing is that bvp5c, bvp4c do is that they minimizes the (2-)norm of the residual by searching ya, yb as control parameters. It must internally the compute the Jacobian of the residual with respect to ya and yb. By chain rule it have to compute the derivative dyb/dya (coupled by the forward ODE and possibly the inverse derivative dya/dyb by the backward ODE. As it handles the internal time finite-difference discretization, it can compute the Jacobian by chain rule.

The infamous error message "Jacobian singular" occurs when it detects non sensitivity of the residual to the control states ya/yb, and the Newton-like method must inverse singular Jacobian, which is not possible numerically.

That's why I think correctly normalization of the statevector and it's derivative are very crucial for bvp5c, bvp4c. The accumulation chain-rule to come to Jacobian will amplify any defect of non-ideal normalization, and ends up having a Jacobian that becomes numerically singular. The sixth order non-linear system make me very nervuous about the feasibility to solve correctly such system.