24 views (last 30 days)
Show older comments
J on 29 Jul 2024 at 8:56
Commented: Torsten on 29 Jul 2024 at 15:23
Accepted Answer: Torsten
- W2.m
Open in MATLAB Online
I tried to solve the bvp using bvp4c procedure. But I have an error "Unable to solve the collocation equations -- a singular Jacobian encountered." How to resolve this issue?
proj()
Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in solution>proj (line 31)
sol= bvp4c(@projfun,@projbc,solinit,options);
function sol= proj
A = 0.5;
Gr = 0.7;
Gc = 0.5;
Kp = 3.0;
beta = 0.5;
Pr = 0.3;
Df = 0.2;
Sc = 0.1;
L0 = 0.5;
Sr = 0.3;
M = 1;
Bi1 = 0.5;
Bi2 = 0.5;
K0 = 0.3;
myLegend1 = {};myLegend2 = {};
rr = [0.3 0.5 0.8];
for i =1:numel(rr)
Re = rr(i);
y0 = [1, 0, 1, 1, 0, 1, 0];
options =bvpset('stats','on','RelTol',1e-4);
x = linspace(0,10,500);
solinit = bvpinit(x,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
disp((sol.y(1,20)))
figure(1)
plot(sol.x,(sol.y(6,:)))
grid on,hold on
myLegend1{i}=['Pr = ',num2str(rr(i))];
xlabel('eta');
ylabel('(thetas-thetaf)/thetas');
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
function dy= projfun(~,y)
dy= zeros(7,1);
dy(1) = y(2);
dy(2) = y(3);
dy(3) = (0.5*A*y(3) - Gr*y(4) - Gc*y(6) + (M + Kp + A)*y(2)) / (1 + (1/beta));
dy(4) = y(5);
dy(5) = (Pr*(0.5*A*y(5)) - Pr*Df*Sc*(0.5*A*y(7) + 2*A*y(6) - L0*y(6))) / (1 - (Pr*Df*Sc*Sr));
dy(6) = y(7);
dy(7) = (Sc*(0.5*A*y(7) + 2*A*dy(6)) - Sr*Pr*(0.5*A*y(5) + 2*A*y(4))) / (1 - (Sr*Df*Pr));
end
function res= projbc(ya,yb)
res= [
ya(2) - (1 + K0*ya(3));
ya(5) + Bi1*(1 - ya(4));
ya(7) + Bi2*(1 - ya(6));
yb(2);
yb(4);
yb(6);
yb(7)];
end
end
1 Comment Show -1 older commentsHide -1 older comments
Show -1 older commentsHide -1 older comments
Torsten on 29 Jul 2024 at 9:11
Direct link to this comment
https://ms-www.mathworks.com/matlabcentral/answers/2141026-how-can-i-solve-the-error-unable-to-solve-the-collocation-equations-a-singular-jacobian-encounte#comment_3222876
How to resolve this issue?
It's not a technical issue. Usually, there is something wrong with the equations or the boundary conditions. So you should compare both of them with the mathematical description of the problem.
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Torsten on 29 Jul 2024 at 13:53
Edited: Torsten on 29 Jul 2024 at 14:53
Open in MATLAB Online
Your ODE system is linear - thus you can determine its general solution having 7 free parameters. But incorporating of your boundary conditions leads to a linear system of equations where the coefficient matrix A is rank-deficient (rank 6 instead of rank 7). Consequently, your system is not solvable.
A = 0.5;
Gr = 0.7;
Gc = 0.5;
Kp = 3.0;
beta = 0.5;
Pr = 0.3;
Df = 0.2;
Sc = 0.1;
L0 = 0.5;
Sr = 0.3;
M = 1;
Bi1 = 0.5;
Bi2 = 0.5;
K0 = 0.3;
Mat = zeros(7);
Mat(1,2) = 1;
Mat(2,3) = 1;
Mat(3,2) = (M + Kp + A) / (1 + 1/beta);
Mat(3,3) = 0.5*A/(1 + 1/beta);
Mat(3,4) = -Gr/(1 + 1/beta);
Mat(3,6) = -Gc/(1 + 1/beta);
Mat(4,5) = 1;
Mat(5,5) = Pr*0.5*A/(1 - Pr*Df*Sc*Sr);
Mat(5,6) = - Pr*Df*Sc*(2*A - L0) / (1 - Pr*Df*Sc*Sr);
Mat(5,7) = - Pr*Df*Sc*0.5*A/ (1 - Pr*Df*Sc*Sr);
Mat(6,7) = 1;
Mat(7,4) = - Sr*Pr*2*A / (1 - Sr*Df*Pr);
Mat(7,5) = - Sr*Pr*0.5*A / (1 - Sr*Df*Pr);
Mat(7,7) = Sc*(0.5*A+2*A) / (1 - Sr*Df*Pr);
syms x1(t) x2(t) x3(t) x4(t) x5(t) x6(t) x7(t)
eqns = [diff(x1);diff(x2);diff(x3);diff(x4);diff(x5);diff(x6);diff(x7)]-Mat*[x1;x2;x3;x4;x5;x6;x7]==[0;0;0;0;0;0;0];
sol = dsolve(eqns,'MaxDegree',4)
sol = struct with fields:
x2: C3*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((... x1: C1 + C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((33648505967... x3: C6*exp(-(t*(865^(1/2) - 1))/24) + C7*exp((t*(865^(1/2) + 1))/24) - C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 23213488... x4: C2*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((... x5: C4*exp((t*(9802324*((3364850596725*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/3))/960855558009... x6: C3*(exp((t*(2976375*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/6) - 9802324*((3364850596725*((... x7: C4*exp((t*(9802324*((3364850596725*((13875*3^(1/2)*342589119285075140699353^(1/2))/1883723499362492536448 - 2321348855894375/1076413428207138592256)^(1/3))/960855558009...
eqn1 = subs(sol.x2,t,0)-(1+K0*subs(sol.x3,t,0))==0;
eqn2 = subs(sol.x5,t,0)+Bi1*(1-subs(sol.x4,t,0))==0;
eqn3 = subs(sol.x7,t,0)+Bi2*(1-subs(sol.x6,t,0))==0;
eqn4 = subs(sol.x2,t,10)==0;
eqn5 = subs(sol.x4,t,10)==0;
eqn6 = subs(sol.x6,t,10)==0;
eqn7 = subs(sol.x7,t,10)==0;
[A,b] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7])
A=
b=
rank(A)
ans = 6
rank([A,b])
ans = 7
1 Comment Show -1 older commentsHide -1 older comments
Show -1 older commentsHide -1 older comments
Torsten on 29 Jul 2024 at 15:23
Direct link to this comment
https://ms-www.mathworks.com/matlabcentral/answers/2141026-how-can-i-solve-the-error-unable-to-solve-the-collocation-equations-a-singular-jacobian-encounte#comment_3223231
I suspect the reason is that you don't have any boundary condition involving y(1).
Sign in to comment.
More Answers (0)
Sign in to answer this question.
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
