bvp4c singular jacobian error

Hello,

I am trying to solve a boundary value problem. But I have faced with " Singular Jacobian Error". Does anybody know how I can solve it?

It doesn't work for Nmf>4 !!!

clc; close all; clear all;
format long e

solinit = bvpinit(x,@xinitfun);
sol= bvp4c(@odefun,@bcfun,solinit);

x= linspace(-1,1,100)  
function xinit= xinitfun(x)
xinit=rand((2*Nmf+1)*6 * 2,1);

function dydx = odefun(x , y)  
global Nmf Re alpha Ra Pr NR

dydx = zeros((2 * Nmf + 1) * 2 * NR,1);
dydx1= zeros((2 * Nmf + 1) * NR,1);

for ii = 1 : 2 * Nmf + 1

      RI = (ii - 1) * NR;
      II = (ii - 1) * NR + (2 * Nmf +1) * NR;
      n = ii - Nmf - 1;
      na = n * alpha;

      teta_f2 = 0;
      NLT1 = 0;
      NLT2 = 0;

      for jj = 1 : 2 * Nmf + 1

          m = jj - Nmf - 1;
          ma = m * alpha;
          nma =(n - m) * alpha;
          teta_f = 0;
          diff_teta_f = 0;

          if abs(n - m) <= Nmf

              if abs(n - m) == 1

                  teta_f = (-1 / 8 * Pr) * sinh(nma .* x) ./ sinh(nma) + (1 / 8 * Pr) * cosh(nma .* x) ./ cosh(nma);
                  diff_teta_f = (nma * sinh(nma .* x)) ./ (8 * Pr * cosh(nma)) - (nma * cosh(nma .* x)) ./ (8 * Pr * sinh(nma));

              elseif (n - m) == 0

                  teta_f = 0.5 .* (1 - x) / Pr;
                  diff_teta_f = 0.5 / Pr;

              end

              NLT1 = NLT1 + (y(2 + RI) + y(2 + II) * i) * ma .* ((y(3 + RI) + y(3 + II) * i) - ma ^ 2 .* (y(1 + RI) + y(1 + II) * i))...
                  - nma .* (y(1 + RI) + y(1 + II) * i) .* ((y(4 + RI) + y(4 + II) * i) - ma ^ 2 .* (y(2 + RI) + y(2 + II) * i));

              NLT2 = NLT2 + nma .* (y(2 + RI) + y(2 + II) * i) .* (teta_f + Pr .* (y(5 + RI) + y(5 + II) * i)) - ma .*...
                  (y(1 + RI) + y(1 + II) * i) .* (diff_teta_f + Pr .*( y(6 + RI) + y(6 + II) * i));

          end

      end

      if abs(n) == 1

          teta_f2 = (-1 / 8 * Pr) * sinh(na .* x) ./ sinh(na) + (1 / 8 * Pr) * cosh(na .* x) ./ cosh(na);

      elseif n == 0

          teta_f2 = 0.5 .* (1 - x)/ Pr;

      end

      dydx1((ii - 1) * NR + 1 : ii * NR,1) = [y(2 + RI) + y(2 + II) * i
          y(3 + RI) + y(3 + II) * i
          y(4 + RI) + y(4 + II) * i
          2 * (na) ^ 2 .* (y(3 + RI) + y(3 + II) * i) - na ^ 4 .* (y(1 + RI) + y(1 + II) * i) + (na * Re .* (1 - x .^ 2) .* ...
          (y(3 + RI) + y(3 + II) * i - na ^ 2 .* (y(1 + RI) + y(1 + II) * i)) - na * Re * (-2) .* (y(1 + RI) + y(1 + II) * i) + na * Ra .* (y(5 + RI) + y(5 + II) * i) + (Ra/Pr) * na .* teta_f2 + NLT1) * i
          y(6 + RI) + y(6 + II) * i
          na ^2 .* (y(5 + RI) + y(5 + II) * i) + (na * Re * (1 - x .^ 2) .* (Pr .* (y(5 + RI) + y(5 + II) * i) + teta_f2 ) + NLT2) * i];

end

dydx=[real(dydx1);imag(dydx1)];

    function res=bcfun(ya,yb)
    global Nmf NR

res=[ya(1 : NR :  (2 * Nmf + 1) * 2 * NR)
     ya(2 : NR :  (2 * Nmf + 1) * 2 * NR)
     ya(5 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(1 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(2 : NR :  (2 * Nmf + 1) * 2 * NR)
     yb(5 : NR :  (2 * Nmf + 1) * 2 * NR)];

Answer 1

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/6369-bvp4c

Use the SingularTerms option.

[Amir K Neghab]

Thanks Walter. I have found my silly errors.

Instead of y(6) I had to write y(6+RI).

But i would like to get to know what "SingularTerms option"is.

Thanks so much for your prompt answer

Amir

Answer 2

[Walter Roberson]

In your revised code, the "x=" line must go before the "solinit=" line.

When you say that "It doesn't work for Nmf>4", is it that you see the message about singular jacobian, or is it that you see something else happen? And for certainty, Nmf must be a positive integer right? And does the code fail when Nmf is exactly equal to 4, or only when Nmf is greater than 4 (that is, 5 or higher) ?

You have a number of global variables whose values we do not know, so we cannot test your code.

[Amir K Neghab]

Dear Walter,

Thanks for responding.

1. "x=" line is before the "solinit=" in my file

2. yes, Nmf is a positive integer.

3. The code gives me singular jacobian error when Nmf is >=4. I mean it only works for Nmf=1 (mode -1,0 and+1) Nmf=2 (Mode -2,-1,...,+1&2) and Nmf=3 (Mode -3,-2,...,+3)

4. The global variables are :

Nmf=1; % Number of Fourier Modes of Mean Flow

alpha=5;

Re=0.05;

Pr=0.71;

Ra=2000;

NR=6;

6. if you need more info I can send you my code to your email address.

Thanks again for your time