6th order differential equation [bvp5c issue]

10 Ansichten (letzte 30 Tage)
Wissem-Eddine KHATLA
Wissem-Eddine KHATLA am 26 Apr. 2022
Bearbeitet: Bruno Luong am 27 Apr. 2022
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.

Antworten (2)

Torsten
Torsten am 26 Apr. 2022
Bearbeitet: Torsten am 26 Apr. 2022
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.
  7 Kommentare
Wissem-Eddine KHATLA
Wissem-Eddine KHATLA am 26 Apr. 2022
@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
Torsten am 26 Apr. 2022
Nevertheless, you should give it a try. My attempt was not very thorough.

Melden Sie sich an, um zu kommentieren.


Wissem-Eddine KHATLA
Wissem-Eddine KHATLA am 26 Apr. 2022
@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.
  2 Kommentare
Torsten
Torsten am 26 Apr. 2022
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
Bruno Luong am 26 Apr. 2022
Bearbeitet: Bruno Luong am 27 Apr. 2022
@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.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by