2

I am trying to solve this differential equation : enter image description here

In order to do so : there is a mathematical part consisting on determining the boundary conditions of this problem that are defined as :

enter image description here

More precisely, the first BC is formulated differently since we need an extra condition in order to "close" the system.

We choose to write : enter image description here. Then, we solve asymptotically the new problem : enter image description here. Therefore, there are 5 roots as a solution to this differential equation : 2 of these roots are growing exponentials that the following condition allows to eliminate :

enter image description here

We deduce two BC written as :

enter image description here enter image description here

In addition to the BC :

enter image description here

And : F(0) = 1

I have tried to implement this problem by using bvp5c of Matlab but I haven't succeeded. Here is my code :

eta0 = 20;
etamesh = [linspace(-eta0,0,100),linspace(0,eta0,100)];
Finit = [1; 0.5; 1; 0; 0];

solinit = bvpinit(etamesh,Finit)
sol = bvp5c(@fun_ode, @bc, solinit);

figure()
plot(sol.x,sol.y(1,:),'linewidth',1.5) % Plotting F
ylim([-2 2]);
hold on
plot(sol.x,sol.y(3,:),'linewidth',1.5) % Plotting F''
plot(sol.x,sol.y(5,:),'linewidth',1.5) % Plotting F''''
grid on

legend("F","F''","F''''",'location','northwest')

function dFdeta = fun_ode(eta,F,region)
  dFdeta=[F(2);F(3);F(4);F(5);(1-F(1))/F(1)^3];
end

% Boundary conditions : 
% F'''(-eta0) = F''''(-eta0) = 0
% F(0) = 1
% [1] = [2] = 0 in eta = eta0


function res = bc(FL,FR)

alpha = cos(3*pi/5) + 1i*sin(3*pi/5);
beta  = cos(7*pi/5) + 1i*sin(7*pi/5);
gamma = abs(alpha);

  res = [FL(4,1) ; FL(5,1);...

         FR(1,1) - 1 ;...
         FR(1,1)-FL(1,2) ; FR(2,1)-FL(2,2) ; FR(3,1)-FL(3,2) ; FR(4,1)-FL(4,2) ;  FR(5,1)-FL(5,2) ;...
         
         (1-(beta+alpha))*FR(3,2) + (-(alpha+beta)+gamma^2)*FR(2,2) ; (-(beta+alpha)+gamma^2)*FR(3,2) + (FR(2,2))*gamma^2];
end

I am mainly struggling with the implementation of the two coupled BC written above and the continuity condition for F(0)...

I have commented my code and I would be grateful if someone give some time to look for potential errors or providing some advices.

Thank you very much.

Wiss
  • 145
  • 6
  • In general this looks sensible, but in the code so far you have not the values at 3 points for the conditions. I vaguely remember that the matlab bvp solvers allow such multi-point conditions, but that will surely involve setting options. // Alternatively you could use a single-shooting approach, set the derivative values at `x=0` as variable, integrate in both directions to evaluate the boundary conditions and return their residuals, put that function into `fsolve`. – Lutz Lehmann Jun 06 '22 at 12:10
  • @LutzLehmann Thank you for your reply : I have never done that before that is why I have tried to implement my code like shared in my post. Do you advise using options ? Are you familiar with the procedure that you are suggesting ? Because it would be easy for me to take inspiration from an existing test case. Thank you – Wiss Jun 06 '22 at 14:29

0 Answers0