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

I've posted yesterday a question about solving a non linear equation : it was not clear so I am reformulating my question.

I am trying to solve a high-order non linear differential equation presented right below with $$F \in [-\eta_{0}, \eta_{0}]$$

$$\left(5F - 4\eta\frac{dF}{d\eta}\right) = 9 \frac{d}{d\eta}\left(F^{3} \frac{d^{5}F}{d\eta^{5}}\right)$$

As boundary conditions, I have the following ones :

$$\left\{ \begin{array}{ll} F' = F'''=0 & \mbox{in } \space \eta = 0 \\ F = F' = 0 & \mbox{in} \space \eta =\pm \space\eta_{0} \end{array} \right.$$

I am trying to apply a shooting method algorithm in order to solve this equation on $$[0, \eta_{0}]$$, I have converted it in a system of 1-st order differential equations as following :

$$f_{i}' = \left[\begin{array}{@{}c@{}} y_{2}=\frac{dF}{d\eta} \\ y_{3} = \frac{d^{2}F}{d\eta^{2}}\\ y_{4} = \frac{d^{3}F}{d\eta^{3}} \\ y_{5} = \frac{d^{4}F}{d\eta^{4}} \\ y_{6} = \frac{d^{5}F}{d\eta^{5}}\\ y_{7} = \frac{d^{6}F}{d\eta^{6}} \end{array} \right]$$ with : $$y_{7} = \left[\left(\frac{5}{9y_{1}^{2}}\right) - \left(\frac{4xy_{2}}{9y_{1}^{3}}\right)-\left(\frac{3y_{2}y_{6}}{y_{1}}\right)\right]$$

According to the boundary conditions we thus have : $$y_{2}(0) = y_{4}(0) = 0$$ and we have to guess the remaining ones : $$y_{1}(0), y_{3}(0),y_{5}(0) \space \mbox{and} \space y_{6}(0)$$

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 value of $$\eta_{0}$$ very "high" in order to match the conditions...

Here are the algorithms that I am using : These are MATLAB scripts so I am sharing it below. The main function is the one called "shoot4nl".

function shoot4nl

% Shooting method for nonlinear 4th-order boundary
% value problem

global XSTART XSTOP H   % Make these params. global.
XSTART = 0; XSTOP = 50; % Range of integration.
H = 0.25;                % Step size.
freq = 1;               % Frequency of printout.
u = [1 -1 1 0];             % Trial values of u(1)
% and u(2), u(3) and u(4).

x = XSTART;
u = newtonRaphson2(@residual,u);
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
printSol(xSol,ySol,freq)

function F = dEqs(x,y)    % Differential equations.
F = zeros(1,6);
F(1) = y(2); F(2) = y(3); F(3) = y(4); F(4) = y(5); F(5) = y(6);
F(6) = ((5/(9*y(1)^2))-(4*x*y(2)/(9*y(1)^3))-(3*y(2)*y(6)/y(1)));

function y = inCond(u)    % Initial conditions; u(1),u(2),u(3) and u(4) are unknowns
y = [u(1) 0 u(2) 0 u(3) u(4)];

function r = residual(u)  % Bounday residuals.
global XSTART XSTOP H
r = zeros(length(u),1);
x = XSTART;
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
lastRow = size(ySol,1);
r(1) = ySol(lastRow,3);
r(2) = ySol(lastRow,4) - 1;


In order to run it we have to define other functions.

runKut5 :

function [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)

% 5th-order Runge-Kutta integration.
% USAGE: [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)
% INPUT:
% dEqs  = handle of function that specifyies the 1-st order differential
% equations : F(x,y) = [dy1/dx dy2/dx dy3/dx ...]

%
% x,y
% xStop = terminal value of x.
% h     = trial value of increment of x.
% eTol  = per-step error tolerance (default = 1.0e-6).
% OUTPUT:
% xSol = x-values at which solution is computed.
% ySol = values of y corresponding to the x-values.

if size(y,1) > 1 ; y = y'; end % y must be row vector

if nargin < 6; eTol = 1.0e-6; end

n = length(y);
A = [0 1/5 3/10 3/5 1 7/8];
B=[    0          0       0          0           0
1/5        0       0          0           0
3/40      9/40     0          0           0
3/10     -9/10    6/5         0           0
-11/54      5/2   -70/27      35/27         0
1631/55296 175/512 575/13824 44275/110592 253/4096];
C = [37/378 0 250/621 125/594 0 512/1771];
D = [2825/27648 0 18575/48384 13525/55296 277/14336 1/4];
% Initialize solution
xSol = zeros(2,1); ySol = zeros(2,n);
xSol(1) = x; ySol(1,:) = y;
stopper = 0; k = 1;
for p = 2:5000
% Compute K’s from Eq. (7.18)
K = zeros(6,n);
K(1,:) = h*feval(dEqs,x,y);
for i = 2:6
BK = zeros(1,n);
for j = 1:i-1
BK = BK + B(i,j)*K(j,:);
end
K(i,:) = h*feval(dEqs, x + A(i)*h, y + BK);
end
% Compute change in y and per-step error from
% Eqs.(7.19) & (7.20)
dy = zeros(1,n); E = zeros(1,n);
for i = 1:6
dy = dy + C(i)*K(i,:);
E = E + (C(i) - D(i))*K(i,:);
end
e = sqrt(sum(E.*E)/n);
% If error within tolerance, accept results and
% check for termination
if e <= eTol
y = y + dy; x = x + h;
k = k + 1;
xSol(k) = x; ySol(k,:) = y;
if stopper == 1;
break
end
end
% Size of next integration step from Eq. (7.24)
if e~=0; hNext = 0.9*h*(eTol/e)^0.2;
else; hNext=h;
end
% Check if next step is the last one (works
% with positive and negative h)
if (h > 0) == (x + hNext >= xStop )
hNext = xStop - x; stopper = 1;
end
h = hNext;
end


newtonRaphson2 :

function root = newtonRaphson2(func,x,tol)

% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
% USAGE: root = newtonRaphson2(func,x,tol)
% INPUT:
% func = handle of function that returns[f1,f2,...,fn].
% x    = starting solution vector [x1,x2,...,xn].
% tol  = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.

if nargin == 2
tol = 1.0e4*eps;
end

if size(x,1) == 1; x = x'; end   % x must be column vector
for i = 1:30
[jac,f0] = jacobian(func,x);
if sqrt(dot(f0,f0)/length(x)) < tol
root = x; return
end
dx = jac\(-f0);
x = x + dx;
if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0)
root = x; return
end
end
error('Too many iterations')
function [jac,f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).
h = 1.0e-4;
n = length(x);
jac = zeros(n);
f0 = feval(func,x);
for i =1:n
temp = x(i);
x(i) = temp + h;
f1 = feval(func,x);
x(i) = temp;
jac(:,i) = (f1 - f0)/h;
end


and printSol :

function printSol(xSol,ySol,freq)

% Prints xSol and ySoln arrays in tabular format.
% USAGE: printSol(xSol,ySol,freq)
% freq = printout frequency (prints every freq-th
%        line of xSol and ySol).

[m,n] = size(ySol);
if freq == 0;freq = m; end
for i = 1:n
end
for i = 1:freq:m
fprintf('%14.4e',xSol(i),ySol(i,:)); fprintf('\n')
end
if i  ~= m; fprintf('%14.4e',xSol(m),ySol(m,:)); end


Unfortunately, when I run it : I have an error stating that one of the matrices is singular : I assume that there is somewhere a division by 0 that occurs...

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


Could someone help me with his mathematical expertise in order to solve this problem ?

• Are we looking for an even solution here, $F(\eta)=F(-\eta)$? Apr 7 at 20:06