unable to identify the error in my program so as to get my plots
조회 수: 1 (최근 30일)
이전 댓글 표시
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
];
댓글 수: 0
채택된 답변
Walter Roberson
2021년 2월 7일
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
Walter Roberson
2021년 2월 11일
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.
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!