numerical methods - Lack of 'convergence' in MATLAB's bvp5c

Question

I am trying to solve the following set of 2nd order non-linear and coupled ODEs:

    0 = ??2 ?″(??) ? ?? ?′(??) + ??2 ??(??)2 [1??(??)], 
    0 = ??2 ??″(??) + ?? ??'(??) ? ?? ??2 ??(??) [??(??)2 + ??(??)2 ? 2], 
    0 = ??2 ??″(??) + ?? ??′(??) ? (1/2) ??(??) [1??(??)]2 ? ?? ??2 ??(??) [??(??)2 + ??(??)2 ? 2].

(I'm sorry EQs do not look beautiful. I'm used to LaTeX syntax only). Well, besides the equations, I have the following Boundary Conditions:

    f'(0) = 0, f'(x → ∞) = 0. I also know that f (x → ∞) = 1 ,
    h(0)  = 0, h(x → ∞)  = 1 ,
    g(0)  = 0, g(x → ∞)  = 1 . 

Moreover, I also expect the first derivatives h'(x), f'(x), and g'(x) to go to zero for some finite value of x and then, just stay zero. That is, I know my solution must reach the asymptotic values and remain constant afterwards. In other words, I know the functions h(x), f(x) and g(x) must 'saturate'.

Using MATLAB, I have tried the following solution:

xmin=1e-3;
xmax=50;
guess = [ 1 1 1 0 0 0];


xmesh = linspace(xmin,xmax,1e5);
solinit = bvpinit(xmesh,guess);%The last vector is my guess.
options = bvpset('RelTol',1e-5,'NMax',5e6); 

sol = bvp5c(@deriv, @bcs, solinit, options);

Sxint = deval(sol,xmesh);
figure
plot(sol.x(1,:),sol.y(1,:),'k-');
hold on 
plot(sol.x(1,:),sol.y(2,:),'m-');
hold on
plot(sol.x(1,:),sol.y(3,:),'k-')
hold off
axis([0 xmax -0.2 1.5])

function dydx = deriv(x,y)
lambda=0.5;
dydx= [ y(4) %The vector y() was keeping the following results: y=(h, f, g, h', f', g')
        y(5)
        y(6)
        (1/x)*y(4) - (y(3)^2)*(1-y(1))
        -(1/x)*y(5) + (lambda)*y(2)*(y(2)^2 + y(3)^2 - 2)
        -(1/x)*y(6) + (1/(2*x^2))*y(3)*((1-y(1))^2) + (lambda)*y(3)*(y(2)^2 + y(3)^2 - 2)];
end

% boundary conditions 

function res = bcs(ya,yb)
res = [ya(1)
       yb(1) - 1
       yb(2) - 1
       ya(3)
       yb(3) - 1
       ya(5)];
end

Well, up to some minor typos I could make while copying and pasting my code, this code gives me a solution that only takes the desired value (of 'infinity') at the boundary. I can use larger and larger values of my xmax and even so, the solution never starts to saturate at its value for infinity.

I tried using a better guess, based on the analytical solution of these equations for a small value of x, but nothing is giving me a good solution. And this is the reason I am asking for some advice here. What do you think? Is MATLAB incapable of solving this because of the 1/x2 in the third ODE?

Thanks!

question from:https://stackoverflow.com/questions/65928938/lack-of-convergence-in-matlabs-bvp5c

Answer 1

At a first glance, for large x you get in f,g a mechanical system with Hamiltonian

H = 0.5·(f'(x)2+g'(x)2) - 0.25·??·(f(x)2+g(x)2-2)2

with a friction that asymptotically falls to zero. The potential energy is a hill with sides going to -oo and a valley at the circular top, like an old, weathered volcano.

The initial conditions say that the movement starts in the center of the indentation. With your current approach, you reach the rim around the indentation in finite time, and thus with non-zero velocity. I would expect solutions to continue from that to fall down along the outside of the potential energy landscape.

You are trying to reach f=g=1 asymptotically. This corresponds to the energy surface H=0. Intuition (which could be wrong) and angular momentum, which has to converge to zero, says that the final approach has to be radially, that is, with f(x)=g(x) and thus

f'(x)2 = ??·(f(x)2-1)2 ==> f'(x) = c·(1-f(x)2),  c=sqrt(??)
f(x)=tanh(c*(x-d)) 

This asymptotic behavior suggests that the boundary conditions at infinity get replaced by conditions at x=T

f'(T)-c*(1-f(T)^2)=0
g'(T)-c*(1-g(T)^2)=0
h'(T)+g(T)*(1-h(T))=0

the last to kill the exponentially increasing term of h''+g2·(1-h)=0.

One also needs some sensible approximation at x=0, disregarding cubic and higher order terms one gets the following:

• h(x)=h_2·x2 follows from the first equation, with h_2 free,
• f(x)=f_0+f_2·x2 leads to 4*f_2 = ??·f_0·(f_02-2),
• g(x)=g_1·x+g_2·x2 leads to g_1·x+4·g_2·x2 = (g_1·x+g_2·x2)/2 so that g_1=g_2=0.

Numerical experiments to follow. Or not. Up to now the approximations at the upper boundary are not sufficient to force convergence. The more sensible results are close to the circle segment from angle 0 to pi/4