Boundary value problems for ODE with bvp4c - singular jacobian error

Hello all,

I'm trying to solve a boundary value problem with 2 coupled differential equations using bvp4c. I'm new to that kind of problem so i did the tutorial by Jacek Kierzenka available on file exchange. When i try to apply it to my case, i have the singular jacobian error but i don't get the problem.

Here are the equations i try to solve :

• u''=-alpha.u'.(1-exp(-x)/(x.g)
• g'=beta.u^2/x^2

With alpha and beta 2 constant parameters and the following set of boundary conditions :

• u=0 for x=0
• u and g -> 1 for x-> infinity

Here is my code :

infinity = 4;
solinit=bvpinit(linspace(0,infinity,8),[0 1 1]);
options = bvpset('stats','on');
sol=bvp4c(@funODE,@funBC,solinit,options);

%------
function dydx = funODE(x,y)

beta=0.0217;
alpha=2;

dydx=[y(2)
      -alpha*y(2)*(1-exp(-x))/(x*y(3))
      beta*y(1)^2/x^2];

%------
function res = funBC(y0,yinf)
res= [y0(1)
      yinf(1) - 1
      yinf(3) - 1];

I tried to change infinity value, number of points in solinit, values of the initial guess and values for alpha and beta. In every case, i have the singular jacobian error, but i don't get if the problem comes from my formulation or from the initial guess.

Does somebody see where the problem can be ?

Answer 1

[Walter Roberson]

Maple says that the system has no solutions. The constraint u(0) = 0 forces u(x) = 0 but that conflicts with u(infinity) = 1

But let us make sure we are discussing the same equations, as your notation is not typical and you are missing a ')'.

The equations I used were

diff(u(x), x, x) = -alpha*(diff(u(x), x))*(1-exp(-x))/(x*g(x))
diff(g(x), x) = beta*u(x)^2/x^2

In particular please check whether it should be (1-exp(-x))/(x*g(x)) (as you coded in your attempt) or should be (1-exp(-x)/(x*g(x))) which is the other possible place to put the missing ')'

If the expression for u is correct then I can establish from the finite boundary condition for u(infinity) that u(x) must be a constant.

[Yann]

Thank you Walter for your so quick answer.

Sorry for the missing ')'. The equations you used are correct.

Indeed, the system has no analytical solutions, that's why I need to solve it numerically.

Do you think I can do it with a different set of boundary conditions ?

This problem comes from a book and the author solves it with this set of boundaries. I checked the equations which seem to be correct so maybe i should have a look on the boundaries.

u(x) should definitively be 0 for x=0 but it might have a different value than 1 at infinity.

[Walter Roberson]

It is not a matter of the system being too complex for analytical solution: it is that the system of equations with those boundary conditions is inconsistent and thus will not have either analytical nor numeric solution.

In order for the system to have an alternative solution, u(infinity) must be positive or negative infinity, I figure, and u(x) would be (c * x + 0) for some constant c, where the 0 is the boundary condition u(0) = 0.

There is _perhaps_ room for there to be a singularity at infinity, if g(x) can be convinced to go to 0 at infinity at the same rate as x approaches infinity (leading to 0 * infinity in the denominator for u''(x) )

[Yann]

Ok so it seems I miss something in the problem.

u(0) should definitively be 0 and g(x) should have an asymptotic behaviour, tending toward 1 at infinity. g(0) should be the minimum of g(x).

Does my question have still a place here since it is more a maths problem than a matlab one ?

[Walter Roberson]

Maple solves that as u(x) = 0, g(x) = 1

For the general case without any boundary conditions imposed yet, Maple shows two possibilities, one of which immediately leads to the simple conditions, and the other in terms of u''' . That second possibility is a ratio in which both the numerator and denominator become 0 at x = 0, but because of the derivative terms multiplying the x, it is _perhaps_ possible that the derivatives go to infinity at a rate that "cancel" the 0's. It would not be easy to solve, if it could be solved at all.

[Yann]

Hi Walter,

I finally solved my problem.

As there are x on the denominators of both equations, trying to solve the system for x=0 leads to division by zero.

So, bvp4c is able to solve the system as written on my first post but with a boundary near but not equal to zero.

To solve it, I just changed the initial guess using 1e-6 as first boundary instead of 0 :

solinit=bvpinit(linspace(1e-6,infinity,8),[0 1 1]);

Answer 2

[eh]

Hi I want to solve system of coupler ODEs in ode45 my bondery condition :: Some conditions are initially And others in end . can you helpe mi to enter bondery condition???

[e-mail]