6th order differential equation [bvp5c issue]
9 views (last 30 days)
Show older comments
Edited: Bruno Luong on 27 Apr 2022
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 :
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);
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);
% We plot the result [F; F'] in a new figure
plotting = figure(4)
% We first construct a vector of 1-st order ODE's
% The last one is the constraint term
function dF = fun_ode(eta,F)
% We built residuals of the method with y(b)-a = 0.
function res = bcfun(ya,yb)
res =[ya([2 4 6]); yb(1:3)];
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,
Edited: Torsten on 26 Apr 2022
You divide by F(1)^3.
For the function F to be well-behaved at eta = 1, you must have that
has a zero of order at least 3 at eta = 1.
I don't know whether this is solvable.
Bruno Luong on 26 Apr 2022
Edited: Bruno Luong on 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.
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!