unable to identify the error in my program so as to get my plots

1 view (last 30 days)
function Maxw
global M la Pr Nb Nt W Le ga A B
infinity =10;
% M=0.5; S=0.2; Ld=0.3; R=0.6; N=0.3; Nb=0.5; Nt=0.5; Pr=0.71; Ki=0.5; Le=1;
% Bi=0.5; Ec=0.5; Ws=0.3;
W= 0.3; M=0.1; la=0.5; Pr=0.71; A=0.1; B=0.1; Nt=0.5; Nb=0.5;
Le=0.1;ga=0.01;
Ldva=[0.5 1.0 2.0 3.0];
for i = 1:4
Ld=Ldva(i);
lines = {'b','r','g','c'};
solinit = bvpinit(linspace(0,infinity),[0 1 1 la 1 0 0]);
sol = bvp4c(@shootode,@shootbc,solinit);
eta = sol.x;
f = sol.y;
% fprintf('Nt = %f\t,skin =%f\t, Nus=%f\t, shr=%f\n', M ,f(3,1),-f(5,1),-f(7,1));
% fprintf('%f\t',f(3,1)+(1/2)*Ws*f(3,1)*f(3,1));
% fprintf('%f\t',-f(5,1));
% fprintf('%f\t',-f(7,1));
% %
figure(1)
set(gca,'Fontsize',12);
plot(eta,f(2,:), lines{i},'linewidth',1);
xlabel('\eta');
ylabel('f ^{\prime}(\eta)');
hold on
end
% hold on
% --------------------------------------------------------------------------
function dydx = shootode(eta,f)
%EX5OfE ODE function for Example 5 of the BVP tutorial.
global M Nb Nt Le Pr W A B la ga
dydx = [
f(2)
f(3)
(1/(1+W*f(3)))*(-f(1)*f(3)-f(2)^2+M*f(2)-la^2)
f(5)
(-(Pr*(f(1)*f(5)+Nb*f(5)*f(7)+Nt*f(5)*f(5)))-A*f(2)-B*f(4))
f(7)
(((Nt/Nb)*(Pr*(f(1)*f(5)+Nb*f(5)*f(7)+Nt*f(5)*f(5)))-A*f(2)-B*f(4))+ga*Le*f(6)+Le*f(1)*f(7))];
% --------------------------------------------------------------------------
function res = shootbc(fa,fb)
%EX5BC Bounfarf conditions for Example 5 of the BVP tutorial.
global la
res = [ fa(1)-0
fa(2)-1
fa(4)-1
fa(6)-1
fb(2)-la
fb(4)-0
fb(6)-0
];

Accepted Answer

Walter Roberson
Walter Roberson on 7 Feb 2021
Maxw;
function Maxw
global M la Pr Nb Nt W Le ga A B
infinity =10;
% M=0.5; S=0.2; Ld=0.3; R=0.6; N=0.3; Nb=0.5; Nt=0.5; Pr=0.71; Ki=0.5; Le=1;
% Bi=0.5; Ec=0.5; Ws=0.3;
W= 0.3; M=0.1; la=0.5; Pr=0.71; A=0.1; B=0.1; Nt=0.5; Nb=0.5;
Le=0.1;ga=0.01;
Ldva=[0.5 1.0 2.0 3.0];
options = bvpset('NMax', 1e6);
for i = 1:4
fprintf('starting iter %d\n', i);
Ld=Ldva(i);
lines = {'b','r','g','c'};
solinit = bvpinit(linspace(0,infinity),[0 1 1 la 1 0 0]);
sol = bvp4c(@shootode,@shootbc,solinit, options);
eta = sol.x;
f = sol.y;
% fprintf('Nt = %f\t,skin =%f\t, Nus=%f\t, shr=%f\n', M ,f(3,1),-f(5,1),-f(7,1));
% fprintf('%f\t',f(3,1)+(1/2)*Ws*f(3,1)*f(3,1));
% fprintf('%f\t',-f(5,1));
% fprintf('%f\t',-f(7,1));
% %
figure(1)
set(gca,'Fontsize',12);
plot(eta,f(2,:), lines{i},'linewidth',1);
xlabel('\eta');
ylabel('f ^{\prime}(\eta)');
hold on
fprintf('done iter %d\n', i);
end
end
% hold on
% --------------------------------------------------------------------------
function dydx = shootode(eta,f)
%EX5OfE ODE function for Example 5 of the BVP tutorial.
global M Nb Nt Le Pr W A B la ga
dydx = [
f(2)
f(3)
(1/(1+W*f(3)))*(-f(1)*f(3)-f(2)^2+M*f(2)-la^2)
f(5)
(-(Pr*(f(1)*f(5)+Nb*f(5)*f(7)+Nt*f(5)*f(5)))-A*f(2)-B*f(4))
f(7)
(((Nt/Nb)*(Pr*(f(1)*f(5)+Nb*f(5)*f(7)+Nt*f(5)*f(5)))-A*f(2)-B*f(4))+ga*Le*f(6)+Le*f(1)*f(7))];
% --------------------------------------------------------------------------
end
function res = shootbc(fa,fb)
%EX5BC Bounfarf conditions for Example 5 of the BVP tutorial.
global la
res = [ fa(1)-0
fa(2)-1
fa(4)-1
fa(6)-1
fb(2)-la
fb(4)-0
fb(6)-0
];
end
On MATLAB Online, the result is:
starting iter 1
Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in q738457>Maxw (line 21)
sol = bvp4c(@shootode,@shootbc,solinit, options);
Error in q738457 (line 1)
Maxw;
Note: When I did not the options configured, it was complaining about the number of mesh points needing to be larger.
  3 Comments
Walter Roberson
Walter Roberson on 11 Feb 2021
Eventually the code ends up doing an LU decomposition on a Jocobian matrix of size 18487, but the rank of the L component turns out to be only 18473 . So this is not a case of the jacobian "slipping" into singular due to simple round-off problems with one row: 14 of the rows become redundant. That is not even 1 in 1000, but it is enough.

Sign in to comment.

More Answers (0)

Categories

Find more on Linear Algebra in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!