Nonlinear differential equation with unknown parameter

Hello,

I am writing this question because I have a problem with a differential equation of order 4 with an unknown parameter noted V :

I have 5 boundary conditions to apply :

I am applying a shooting method : I have found a recent discussion (Solving a 6th order non-linear differential equation in Matlab) and I have tried to implement the method for my problem since I found in the documentation that the solver "bvp4c" can take into consideration an unknown parameter :

V = 0.5;
y_init = [1;0;2;0]; % conditions initiales


psi = linspace(0,200,100);
solinit = bvpinit(eta,y_init,V);
sol = bvp4c(@core, @boundary, solinit);
psi = sol.x;
F = sol.y(1,:);

% Affichage du resultat

figure()
plot(eta,F,'--','linewidth',2,'Color','r');

% Parametre ?

fprintf('Unknown parameter is approximately %7.3f.\n',sol.parameters)

function dF = core(psi,F,V)
dF=[F(2:4);
    ((psi*V*F(2))-V*F(1)-(3*F(2)*F(4)*F(1)^2))/F(1)^3];
end

function res = boundary(ya,yb,V)

res = [ya(1)-1; ya(2); ya(4); yb(2)-1; yb(3)];
end

Unfortunately, I donit obtain the good result because i know that the value of V should be 0.88 : also the plotting obtained changes a lot if I change my interval...

If possible I need help.

Have a nice day.

[Torsten]

I am applying a shooting method

No, you are applying a finite-difference method.

[Karl Zero]

@Torsten I have modified the code of the post mentionned to apply it on my own different problem. I have surely made a mistake but where ?

[Torsten]

No, you didn't make a mistake, I guess.

But in the discussion of (Solving a 6th order non-linear differential equation in Matlab), we didn't yet succeed to get a stable solution, and you also have an ODE of degree 4 which makes things quite difficult.

Answer 1

[Torsten]

Seems your problem is easier than the one discussed under (Solving a 6th order non-linear differential equation in Matlab)

format long
bc0 = [1.0]   %guess for H''(0)
V0 = [2.0];   % Guess for V
p = [bc0,V0];
%p=[0.6323   0.7456]  %eta=10
%p=[0.6369   0.7566]  %eta=15
%p=[0.6369   0.7566]  %eta=20
%p=[0.6562   0.8037]  %eta=25
%p=[0.6608   0.8150]  %eta=30
%p=[0.6641   0.8231]  %eta=40
%p=[0.6643   0.8236]  %eta=50
%p=[0.6637   0.8222]  %eta=60
%p=[0.6627   0.8198]  %eta=80
%p=[0.6705   0.8392]  %eta=120
%p=[0.6688   0.8348]  %eta=160
%p=[0.6828   0.8701]  %eta=200
%p=[0.689030205524681   0.886117860845706] %eta=250
%p=[0.689018898033062   0.886088749625027] %eta=275
%p=[0.689014752596432   0.886078071923821] %eta=300
%p=[0.689013817771355   0.886075659784428] %eta=325
%p=[0.689012647439617   0.886072642654987] %eta=350
%p=[0.689012439215454   0.886072105009582] %eta=375
%p=[0.689012345124427   0.886071857377161] %eta=400
%p=[0.689012324301630   0.886071797757327] %eta=450
%p=[0.689012318896960   0.886071781172971] %eta=500
p=[0.689012317389282   0.886071776014365] %eta=550
etaspan = [0 600];
p = fsolve(@(bc)fun(bc,etaspan),p) % Adjusting the initial condition for H''(0) at eta = 0
bc_total = [1;0;p(1);0];
V = p(2);
[eta,H] = ode45(@(eta,H)fun_ode(eta,H,V),etaspan,bc_total);
plot(eta,H(:,1))

% Creating residuals for the shooting method
function res = fun(p,etaspan)
  bc_total = [1;0;p(1);0];
  V = p(2);
  etaspan_ode = [etaspan(1),etaspan(2)];
  [eta,H] = ode45(@(eta,H)fun_ode(eta,H,V),etaspan_ode,bc_total);
  res(1) = H(end,2)-1.0;  % Deviation of H' from 1 at eta = Inf
  res(2) = H(end,3);      % Deviation of H'' from 0 at eta = Inf
end
% System of 1-st order ODE
function dHdeta = fun_ode(eta,H,V)
  dHdeta=[H(2);H(3);H(4);(-V*H(1)+V*eta*H(2)-3*H(1)^2*H(2)*H(4))/H(1)^3];
end

[Karl Zero]

@Torsten Tank you : this works very well. I would like to know why you have introduced this p vector ? I didn't get what was particularly wrong in my previous code ?

[Torsten]

I didn't get what was particularly wrong in my previous code ?

Nothing. But you said it didn't work as you wanted.

Maybe starting with the solution from my code will also make your code give results as expected.

p(1) = H''(0) and p(2) = V in my code.

For these problems, everything depends on good initial guesses for the unknowns.

[Karl Zero]

@Torsten I have a question : I have tried to keep my original code but with doing the same idea that you suggested. I think it's the same but somehow I have NOT been able to converge to the desired solution. Here is the code :

V0 = [1.0]; % V.
H_0 = [2.0]; %  H''(0).

ca = [H_0,V0]; 

%ca = [-216.812,213.037]; % eta = 10 
%ca = [-5502.389,-530316.562]; % eta = 20 
%ca = [0.191,-0.487]; % eta = 30 
%ca = [-0.002,0.0]; % eta = 50 
ca = [0.664,0.822]; % eta = 60 


y_init = [1.0;0;ca(1);0]; % conditions initiales

V = ca(2);

psi = linspace(0,70);
solinit = bvpinit(psi,y_init,V);
sol = bvp4c(@core, @boundary, solinit);
psi = sol.x;
F = sol.y(1,:);

% Plot

figure()
plot(eta,F,'--','linewidth',2,'Color','r');

% V ? 

fprintf('The value of F^{2}(0) is approximately %7.3f.\n',sol.y(3,1))
fprintf('Unknown parameter is approximately %7.3f.\n',sol.parameters)

function dF = core(psi,F,V)

dF=[F(2);F(3);F(4);(-V*F(1)+V*psi*F(2)-3*F(1)^2*F(2)*F(4))/F(1)^3];

end

function res = boundary(ya,yb,V)

res = [ya(1)-1; ya(2); ya(4); yb(2)-1; yb(3)];
end

I really don't understand why : it's very strange for me. I start with the same guess as yours but my methode diverge... CAn you help please ?

[Torsten]

My code says:

When I used psi = linspace(0,10), the solution with my code was H''(0) = 0.6323, V= 0.7456,

when I used psi = linspace(0,20), the solution with my code was H''(0) = 0.6369, V = 0.7566

and so on.

So you can try setting e.g.

psi = linspace(0,10,100);
y_init = [1.0;0;0.6323;0]; % conditions initiales
V = 0.7456;

and see if your code reproduces these values for H'' and V.

[Karl Zero]

@Torsten But the solution at psi = linspace(0,10) : [0.6323; 0.7456] was obtained after the calculation with the guessing : [2.0;1.0] ? Because, for my code when I start with [2.0; 1.0] I obtain [-216.812,213.037] and after eta = 60 my code starts diverging... I don't know why

[Torsten]

But the solution at psi = linspace(0,10) : [0.6323; 0.7456] was obtained after the calculation with the guessing : [2.0;1.0] ?

Yes. But you have the code - you can check it by yourself.

What do you get if you start with [0.6323; 0.7456] instead of [2.0;1.0] for psi = 10 ?

[Karl Zero]

@Torsten I think I'm not explicit : Your code works very well but I try to understand why my previous code don't work.

For example : I don't get the same result as yours with a starting : of [0.6323; 0.7456] for psi = 10. I obtain very anormal values (-22257; -3455) : I don't understand the reason

THe code I'm talking about :

V0 = [1.0]; % V.
H_0 = [2.0]; %  H''(0).
ca = [H_0,V0]; 

%ca = [-216.812,213.037]; % eta = 10 
%ca = [-5502.389,-530316.562]; % eta = 20 
%ca = [0.191,-0.487]; % eta = 30 
%ca = [-0.002,0.0]; % eta = 50 
ca = [0.664,0.822]; % eta = 60 

y_init = [1.0;0;ca(1);0]; % conditions initiales
V = ca(2);
psi = linspace(0,70);
solinit = bvpinit(psi,y_init,V);
sol = bvp4c(@core, @boundary, solinit);
psi = sol.x;
F = sol.y(1,:);

% Plot
figure()
plot(eta,F,'--','linewidth',2,'Color','r');

% V ? 
fprintf('The value of F^{2}(0) is approximately %7.3f.\n',sol.y(3,1))
fprintf('Unknown parameter is approximately %7.3f.\n',sol.parameters)

function dF = core(psi,F,V)
dF=[F(2);F(3);F(4);(-V*F(1)+V*psi*F(2)-3*F(1)^2*F(2)*F(4))/F(1)^3];
end

function res = boundary(ya,yb,V)
res = [ya(1)-1; ya(2); ya(4); yb(2)-1; yb(3)];
end

[Bruno Luong]

@Karl Zero Torsen does not run MATLAB and bvp4c (he uses octave), he cannot answer on problem with your code.

Answer 2

[Bruno Luong]

Solving up to PsiMax = 1000; which returns V = 0.823. Still not recommended to set PsiMax too large, just to show a more robust code by iterating on the psi span (avoid non linearity) and rescaling so that the problem is better conditioned:

function tata
clear
close all

V = 0.5;
Finit = [1;0;2;0]; % conditions initiales

PsiFinal = 1000;
nbiter = round(10+4*log(PsiFinal));
PsiFinaltab = logspace(log10(1),log10(PsiFinal),nbiter);
N = length(Finit);
for k=1:nbiter
    fprintf('iteration %d/%d\n', k, nbiter);
    psi = linspace(0,PsiFinaltab(k),100);
    psiscale = max(PsiFinaltab(k),1);
    xi = psi/psiscale;
    Ginit = Finit .* psiscale.^(0:N-1)';
    solinit = bvpinit(xi,Ginit,V);
    try
        sol = bvp5c(@(xi,G,V) fun_ode(xi,G,V,psiscale), ...
                    @(ya,yb,V) boundary(ya,yb,V,psiscale), solinit);
    catch
        continue
    end
    Finit = sol.y(:,1) ./ psiscale.^(0:N-1)';
    V = sol.parameters;
end
xi = sol.x;
psi = psiscale*xi;
F = sol.y(1,:);

V = sol.parameters;

fprintf('Unknown parameter is approximately %7.3f.\n', V);

% Affichage du resultat
figure()
plot(psi,F,'--','linewidth',2,'Color','r');

end

function dG = fun_ode(xi,G,V,psiscale)
psi = xi*psiscale;
N = length(G);
F = G ./ (psiscale).^(0:N-1)';
dF=[F(2:4);
    psi*V*F(2)/F(1)^3-V/F(1)^2-3*F(2)*F(4)/F(1)];
dG = dF .* (psiscale).^(1:N)';
end

function res = boundary(ya,yb,V,psiscale)
N = length(ya);
ya = ya ./ (psiscale).^(0:N-1)';
yb = yb ./ (psiscale).^(0:N-1)';
res = [ya([1,2,4]); yb(2:3)]-[1;0;0;1;0];
w = [1; 1; 1; 1; 1];
res = res .* w;
end

[Torsten]

@Bruno Luong

Slope of function at PsiFinal should be 1.

Maybe you forgot to rescale your result for F ?

[Bruno Luong]

@Torsten I do rescale.

yb = yb ./ (psiscale).^(0:N-1)'

[Torsten]

@Bruno Luong

From the plot, it seems that the slope at Infinity is only about 5/1000 instead of 1.

[Bruno Luong]

The solution is not find with precision or does not exist, not because of the set-up.