Error using bvp4c Unable to solve the collocation equations -- a singular Jacobian encountered. Error in MHDI_ori (line 54) sol = bvp4c (@OdeBVP, @OdeBC, solinit, options
2 views (last 30 days)
Show older comments
MHDI_ori()
function MHDI_ori
format long g
%Define all parameters
global R M S Pr a C1 P1 K1 C2 K2 P2 C3 P3 K3 C4 K4 P4 H1 H2 H3 H4 Kn
global D1 D2 B1 B2 B3
%Boundary layer thickness & stepsize
etaMin = 0;
etaMax1 = 15;
etaMax2 = 15; %15, 10
stepsize1 = etaMax1;
stepsize2 = etaMax2;
%Input for the parameters
R= 0.1; %Stretching/Shrinking parameter
M=0.2; %Magnetic field parameter (0,0.1,0.2)
%S=1; %Suction parameter
Pr=6.2; %Prandtl number
a=0.1; %phi1- 1st nanoparticle concentration
%%%%%%%%%%%%%%%%%%%%%% 1st nanoparticle properties (Cu) %%%%%%%%%%%%%%%
C1=385; %C=specific heat
P1=8933; %P=rho=density
K1=400; %K=thermal conductivity
%%%%%%%%%%%%%%%%%%%%%% 2nd nanoparticle properties (Al2O3) %%%%%%%%%%%%
C2=765;
P2=3970;
K2=40;
%%%%%%%%%%%%%%%%%%%%%% 3rd nanoparticle properties (TiO2) %%%%%%%%%%%%%%%
C3=686.2; %C=specific heat
P3=4250; %P=rho=density
K3=8.9538; %K=thermal conductivity
%%%%%%%%%%%%%%%%%%%%%% base fluid %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C4=4179;
P4=997.1;
K4=0.613;
%%%%%%%%%%%%%%%%%%%%%%%%%% multiplier %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H1=P1*C1; %pho*cp nanoparticle 1
H2=P2*C2; %pho*cp nanoparticle 2
H3=P3*C3; %pho*cp nanoparticle 3
H4=P4*C4; %pho*cp nanoparticle 4
B1=C1/C4;
B2=C2/C4;
B3=C3/C4;
Kn=((K1+2*K4)-2*a*(K4-K1))/((K1+2*K4)+a*(K4-K1)); %thermal conductivity nano
D1=1-a;
D2=D1^(2.5);
%%%%%%%%%%%%%%%%%%%%%% first solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10); solinit = bvpinit (linspace (etaMin, etaMax1, stepsize1), @OdeInit1);
sol = bvp4c (@OdeBVP, @OdeBC, solinit, options);
eta = linspace (etaMin, etaMax1, stepsize1);
y = deval (sol, eta);
figure(1) %velocity profile
plot(sol.x,sol.y(2,:),'-')
xlabel('\eta')
ylabel('f\prime(\eta)')
hold on
figure(2) %temperature profile
plot(sol.x,sol.y(4,:),'-')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
%Saving the output in txt file for first solution
descris = [sol.x; sol.y];
save 'ainnaz_upper_ori.txt' descris -ascii
%Displaying the output for first solution
fprintf('\nFirst solution:\n');
fprintf('f"(0) = %7.9f\n',y(3)); %reduced skin friction
fprintf('-theta`(0) = %7.9f\n',-y(5)); %reduced local Nusselt number
fprintf('Cfx = %7.9f\n',(1/D1)*y(3));%skin friction
fprintf('Nux = %7.9f\n',-Kn*y(5) );%local Nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%% second solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options = bvpset('stats','off','RelTol',1e-10);
solinit = bvpinit (linspace (etaMin, etaMax2, stepsize2), @OdeInit2);
sol = bvp4c (@OdeBVP, @OdeBC, solinit, options);
eta = linspace (etaMin, etaMax2, stepsize2);
y = deval (sol, eta);
figure(1)
plot(sol.x,sol.y(2,:),'--')
xlabel('\eta')
ylabel('f`(\eta)')
hold on
figure(2)
plot(sol.x,sol.y(4,:),'--')
xlabel('\eta')
ylabel('\theta(\eta)')
hold on
%Saving the output in txt file for first solution
descris = [sol.x; sol.y];
save 'ainnaz_lower_ori.txt' descris -ascii
%Displaying the output for second solution
fprintf('\nSecond solution:\n');
fprintf('f"(0) = %7.9f\n', y(3)); %reduced skin friction
fprintf('-theta`(0) = %7.9f\n',-y(5)); %reduced local Nusselt number
fprintf('Cfx = %7.9f\n', (1/D1)*(y(3))); %skin friction
fprintf('Nux = %7.9f\n', -Kn*y(5)); %local Nusselt number
fprintf('\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%Define the ODE function
function ff = OdeBVP (x, y);%, M, Pr, Kn, S, R, a, D1, D2, B1, H1, H4)
global M Pr Kn S R a D1 D2 B1 H1 H4
ff = [y(2)
y(3)
D2*(D1+a*B1)*(-(y(1)*y(3))+y(2)*y(2)-1-M+M*y(2))
y(5)
-y(1)*y(5)*Pr*(D1+a*H1/H4)/(Kn)
] ;
end
%Define the boundary condition
function res = OdeBC (ya, yb);%, S, R)
global S R
res = [ya(1)
yb(2)-R
ya(4)-1
yb(2)-1
yb(4)
] ;
end
%Setting the initial guess for first solution
function v = OdeInit1 (x);%, S, R)
global S R
v = [0
-1.25
0
0
0
];
end
%Setting the initial guess for second solution
function v1= OdeInit2 (x)
v1 = [exp(-x)
exp(-x)
exp(-x)
exp(-x)
-exp(-x)];
end
0 Comments
Accepted Answer
Torsten
on 1 Nov 2023
You say you want to have yb(2) - R = 0 and yb(2) - 1 = 0 at the same time as boundary conditions. This is not possible.
More Answers (0)
See Also
Categories
Find more on Equation Solving 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!