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
];

채택된 답변

Walter Roberson
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
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개)

카테고리

Help CenterFile Exchange에서 Linear Algebra에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by