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

25 Ansichten (letzte 30 Tage)
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 Kommentare
Bruno Luong
Bruno Luong am 7 Apr. 2022
Bearbeitet: Bruno Luong am 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
Wissem-Eddine KHATLA am 7 Apr. 2022
@Bruno Luong Thanks for your reply. @Torsten provided a correction to my first attempt by using various in-built functions but it seems like I got an error stating that the "fsolve" function is not correct...

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Torsten
Torsten am 6 Apr. 2022
Bearbeitet: Torsten am 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 Kommentare
Torsten
Torsten am 9 Apr. 2022
Bearbeitet: Torsten am 9 Apr. 2022
I suggest you start with the solution obtained from the shooting method as initial guess.
I'm curious whether bvp4c will reproduce the solution.
Wissem-Eddine KHATLA
Wissem-Eddine KHATLA am 10 Apr. 2022
Bearbeitet: Wissem-Eddine KHATLA am 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,

Melden Sie sich an, um zu kommentieren.


Bruno Luong
Bruno Luong am 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 Kommentare
Torsten
Torsten am 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
Wissem-Eddine KHATLA am 20 Apr. 2022
@Torsten F''(-eta0) doesnt equal 1.35. For instance, here is my plot :
Compared to the one of the article :

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