bvp4c singular jacobian error

6 次查看（过去 30 天）
Amir K Neghab
Amir K Neghab 2011-11-17
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)];

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

采纳的回答

Walter Roberson
Walter Roberson 2011-11-17
  1 个评论
Amir K Neghab
Amir K Neghab 2011-11-18
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

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2011-11-30
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.
  1 个评论
Amir K Neghab
Amir K Neghab 2011-12-1
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Converters (High Power) 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by