Boundary value problems for ODE with bvp4c - singular jacobian error

7 次查看（过去 30 天）
Yann
Yann 2012-1-27
回答： eh 2013-11-29
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 ?

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

采纳的回答

Walter Roberson
Walter Roberson 2012-1-27
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.
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2012-1-31
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
Yann 2012-2-2
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]);

请先登录，再进行评论。

更多回答（1 个）

eh
eh 2013-11-29
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???
haph8712@yahoo.com

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by