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

27 次查看（过去 30 天）

显示更早的评论

Wissem-Eddine KHATLA 2022-4-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab

评论： Wissem-Eddine KHATLA 2022-4-20

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Bruno Luong 2022-4-7

编辑：Bruno Luong 2022-4-7

在 MATLAB Online 中打开

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 2022-4-7

@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...

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Torsten 2022-4-6

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab#answer_936070

编辑：Torsten 2022-4-6

在 MATLAB Online 中打开

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 个评论
显示 28更早的评论隐藏 28更早的评论

Torsten 2022-4-8

编辑：Torsten 2022-4-8

Why did you choose to compute 4 residuals in the function ?

As I already tried to explain to you, I start the ode-integrator with the two conditions set at eta = -eta0 (F = F' = 0) and estimates for the remaining four initial conditions for F'',F''',F'''' and F'''''. You defined four other conditions (two at 0, two at eta = eta0 for F) that should hold. So we have four values that we have to nullify, namely F' = F''' = 0 at eta=0 and F = F' = 0 at eta=eta0. These are the four equations that fsolve have to solve. In the res vector, the errors are returned to fsolve and it is prompted to supply better guesses for F''(-eta0),...,F'''''(-eta0) in order to nullify these errors.

I understand the aim of the line :

bc = fsolve(@(bc)fun(bc,etaspan),bc0);

However, I am not familiar with this syntax (anonymous function) : why did you choose to write it this way ?

"etaspan" is needed in "fun" to define the integration interval for the ode integrator. The line

bc = fsolve(@(bc)fun(bc,etaspan),bc0);

has the effect that "etaspan" is passed to "fun" as an additional parameter that can be used in the function. If you don't like this, you will have to change "etaspan" both in the calling function and in "fun" if you make changes to eta0.

I have strong hope that "bvp4c" can solve your problem more stable than the ad-hoc shooting method I supplied. Did you already test the code I gave you in MATLAB ?

Torsten 2022-4-8

编辑：Torsten 2022-4-8

Since I don't have MATLAB available, I can't test.

But this setting (eta0 = 30, F(eta0)=1e-7) is already very difficult for the solvers. Start with less aggressive values (eta0 = 10, F(eta0) = 1e-3 or something). If you succeed, use the solution as initial value for making eta0 larger and F(eta0) smaller.

But also under OCTAVE, the calculations were quite unstable and small changes in eta0 and F(eta0) led to big changes in the height at eta=0. Although you got the error message of a singular Jacobian, you should try bvp4c or bvp5c with the code I gave to you. Solving boundary value problems with finite differencing is in general much more stable than shooting.

Also solving on a finite interval (e.g. [-1:1] by the coordinate transformation eta~ = tanh(eta)) could stabilize the computations.

Somehow I have the feeling that a condition is missing in the problem formulation. If you start integration at x=0, you have three conditions (because I think the solution is symmetric to eta=0), namely F'(0)=F'''(0)=F'''''(0)=0. Additionally, you have the two conditions at Infinity: F(infinity) = F'(Infinity) = 0. But what is the necessary 6th condition ?

Torsten 2022-4-9

编辑：Torsten 2022-4-9

在 MATLAB Online 中打开

Do you think it's the good numerical reformulation for this particular script ?

Are you aware that this is a totally different problem ? The results of the old formulation gave a value in the order of 1e6 in eta = 0. Now you prescribe F = 1 in eta = 0. Do you think that both formulations have anything in common except for the underlying differential equation ?

This should be the correct code for the reformulated problem, but it doesn't converge:

bc0 = [-2.5004e+05  -2.2123e+06   5.6616e+05]
etaspan = [0 10];
bc = fsolve(@(bc)fun(bc,etaspan),bc0); % Adjusting the initial conditions for F''',F'''',F''''' at eta = eta0
bc_total = [1;0;bc(1);0;bc(2);bc(3)];
[eta,F] = ode45(@fun_ode,etaspan,bc_total);
plot(eta,F(:,1))
% Creating residuals for the shooting method
function res = fun(bc,etaspan)
  bc_total = [1;0;bc(1);0;bc(2);bc(3)];
  etaspan_ode = [etaspan(1),etaspan(2)];
  [eta,F] = ode45(@fun_ode,etaspan_ode,bc_total);
  res(1) = F(end,1)-1e-2;  % Deviation of F from 0 at eta = eta0
  res(2) = F(end,2);      % Deviation of F' from 0 at eta = eta0
  res(3) = F(end,3);      % Deviation of F'' from 0 at eta = eta0
end
% System of 1-st order ODE
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

Torsten 2022-4-9

编辑：Torsten 2022-4-9

在 MATLAB Online 中打开

The solution value at eta = 0 is very important so I need a robust way to determine it.

I don't want to discourage you, but experimenting with the old code showed that the value at eta=0 is strongly dependent on the eta0 chosen.

This works under OCTAVE:

bc0=[   3.8309 0.4494 -0.1066];  
etaspan = [0 12];
bc = fsolve(@(bc)fun(bc,etaspan),bc0); % Adjusting the initial conditions for F''',F'''',F''''' at eta = eta0
bc_total = [bc(1);0;bc(2);0;bc(3);0];
options = odeset('InitialStep',1e-9);
[eta,F] = ode45(@fun_ode,etaspan,bc_total,options);
bc
plot(eta,F(:,1))
end
% Creating residuals for the shooting method
function res = fun(bc,etaspan)
bc_total = [bc(1);0;bc(2);0;bc(3);0];
etaspan_ode = [etaspan(1),etaspan(2)];
options = odeset('InitialStep',1e-9);
[eta,F] = ode45(@fun_ode,etaspan_ode,bc_total,options);
res(1) = F(end,1)-1e-5;  % Deviation of F from 0 at eta = eta0
res(2) = F(end,2);      % Deviation of F' from 0 at eta = eta0
res(3) = F(end,3);      % Deviation of F'' from 0 at eta = eta0
end
% System of 1-st order ODE
function dFdeta = fun_ode(eta,F) 
%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

Torsten 2022-4-9

编辑：Torsten 2022-4-9

https://de.mathworks.com/matlabcentral/answers/120468-why-do-i-receive-an-error-about-a-singular-jacobian-when-i-use-the-bvp4c-function-within-matlab

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 2022-4-10

编辑：Wissem-Eddine KHATLA 2022-4-10

在 MATLAB Online 中打开

@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,

请先登录，再进行评论。

Bruno Luong 2022-4-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1689880-solving-a-6th-order-non-linear-differential-equation-in-matlab#answer_938920

在 MATLAB Online 中打开

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 个评论
显示 77更早的评论隐藏 77更早的评论

Bruno Luong 2022-4-10

在 MATLAB Online 中打开

"However, your last plotting provides the good shape (good values ?). "

I'm not sure as you what mean "good shape" and "good value". I don't know the analytic solution, and I don't know how is other 5 derivative values are correctly matched. Again Imposing condition on fifth derivative is just insane. The derivative is very instable operator, take it five/six time and doing calculation on it is hopeless IMO.

"There is a so strong dependance to the choice of "yinit"...

Because of non linearity, I already tell you that.

When more few algorithms tells that the Jacobian is singular, it just indicates that the problem is so illposed and the dependance of the final state (or at least some subspace of it) to the initial state is numerically zero. It migh be due that your problem is badly scale (for large eta0 >> 1 as I state above) or in the intrincic nature of the problem you are trying to solve.

Also, could you please tell me why did you write the residuals this way ?

res =[ya([2 4 6]); yb(1:3)];

From the bc you are telling us, if you miss my comment in the first code here again

res =[ya(2); % F'(0)=0
    ya(4); % F'''(0)=0
    ya(6); % F'''''(0)=0
    yb(1:3)]; % F(eta0)=0, F'(eta0)=0, F''(eta0)=0,

Bruno Luong 2022-4-15

编辑：Bruno Luong 2022-4-15

在 MATLAB Online 中打开

This code doesn't crash with your experimental data, however there is a warning. In fact your system id 6th order (fake 7th order) and you want to impose 6bc + 1 integral constraint, it cannot meet all in general.

clear
close all
%% Solving the problem of elastic droplet's spreading in an elastic regime in the [0;eta0] interval 
% Solving method using bvp4c solver
% Boundary conditions : F'(0) = F'''(0)= F'''''(0) = 0 and F(eta0) = F'(eta0) = F''(eta0) = 0
% Constraint on the volume
eta0 = 50;
A = 3.7926e-7;
Q = 0.4e-9;
%Finit = [1; 0; 0.5079; 0; 0.1194; 0; 1]; % Initial conditions for the bvpinit solver.
Finit = [-1.0513e+04;0;158.4915;0;-3.6654;0;-1.0513e+04];
nbiter = 20; % 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(log10(1),log10(eta0), nbiter);
for k=1:nbiter
    fprintf('iteration %d/%d\n', k, nbiter);
    etascale = A^(1/6);
    eta = linspace(0,eta0tab(k),100);
    xi = eta/etascale;
    Ginit = Finit .* etascale.^[0:5 0]';
    solinit = bvpinit(xi,Ginit);
    try
        sol = bvp5c(@(xi, G) fun_ode(xi,G,A,Q,etascale), ...
                    @(ya,yb) bcfun(ya,yb,A,Q,etascale), ...
                    solinit);
    catch
%         if k==nbiter
%             fprintf('Fail at the last iteration\n');
%             return
%         end
        continue
    end
    Finit = sol.y(:,1) ./ etascale.^[0:5 0]';
end
xi = sol.x;
eta = etascale*xi;
F = sol.y(1,:);
Fp = sol.y(2,:) / etascale;
% We plot the result [F; F'] in a new figure
figure()
plot(eta,-F,'-','LineWidth',2);
grid on
hold on
plot(eta,Fp,'LineWidth',2)
legend('F',"F'",'Location','northeast')
xlabel('\eta','FontSize',13);
ylabel("F(\eta),F'(\eta)",'fontsize',13)
% We first construct a vector of 1-st order ODE's
% The last one is the constraint term
function dG = fun_ode(xi,G,A,Q,etascale)
eta = xi*etascale;
F = G ./ (etascale).^[0:5 0]';
dF=[F(2:6);
    ((1/(9*A))*(5*F(1)-4*eta*F(2))-3*F(1)^2*F(2)*F(6))/F(1)^3;
    F(1)/A];
dG = dF .* (etascale).^[1:6 1]';
end
% We built residuals of the method with y(b)-a = 0.
function res = bcfun(ya,yb,A,Q,etascale)
ya = ya ./ (etascale).^[0:5 0]';
yb = yb ./ (etascale).^[0:5 0]';
res =[ya([2 4 6]); yb(1:3); yb(7)-Q/A];
end