Coupled Mass and energy equations using bvp4c

5 views (last 30 days)
clc,clear,close
global Pem Peh Da b U yw yf
Pem = 2000; Peh = 500;
Da = 1e8; b = 0.015;
U = 1; yw = 0.9; yf = 0.0503;
z = linspace(0,1);
%guess = 4*ones(1,2);
solinit = bvpinit(z,[1,0,1,0]);
sol = bvp4c(@myfun,@mybc,solinit);
Error using bvp4c (line 197)
Unable to solve the collocation equations -- a singular Jacobian encountered.
function F = myfun(~,y)
global Pem Peh Da b U yw
y1 = y(1);
y2 = y(2);
y3 = y(3);
y4 = y(4);
F(1) = y2;
F(2) = Peh*(y2 + b*Da*y3*exp(-1/y1) - U*(yw-y1));
F(3) = y4;
F(4) = Pem*(y4 + Da*y3*exp(-1/y1));
F = F(:);
end
function f = mybc(ya,yb)
global Pem Peh yf
f(1) = ya(2)+Peh*(yf-ya(1));
f(2) = yb(2);
f(3) = ya(4) + Pem*(1-ya(3));
f(4) = yb(4);
f=f(:);
end
function g = guess(z) % initial guess for y and y'
g = [0.5
0.5
0.5
0.5];
end
I wrote code like this but I'm getting error as
Error using bvp4c (line 248)
Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in Q2 (line 11)
sol = bvp4c(@myfun,@mybc,solinit)
can someone help in this regard
thanks in advance

Accepted Answer

Star Strider
Star Strider on 12 Dec 2021
Ths usual solution for that is to simply increase the initial mesh size, and here that is:
z = linspace(0,1,7500);
It then runs without error. (I used yyaxis because is more than an order of magnitude greater than the others, and those details get lost plotting them on the same axes.)
Also, please never use global variables! Pass them as extra parameters, as I did here.
Increasing the size of ‘z’, correcting the global variable problem, and plotting with yyaxis are the only changes I made to the posted code.
% global Pem Peh Da b U yw yf
Pem = 2000; Peh = 500;
Da = 1e8; b = 0.015;
U = 1; yw = 0.9; yf = 0.0503;
z = linspace(0,1,7500);
%guess = 4*ones(1,2);
solinit = bvpinit(z,[1,0,1,0]);
sol = bvp4c(@(t,y)myfun(t, y, Pem,Peh, Da, b, U, yw, yf),@(ya,yb)mybc(ya,yb, Pem, Peh, yf),solinit)
Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 7500 points and the solution are available in the output argument.
The maximum residual is 2.21838, while requested accuracy is 0.001.
sol = struct with fields:
solver: 'bvp4c' x: [0 1.3335e-04 2.6670e-04 4.0005e-04 5.3340e-04 6.6676e-04 8.0011e-04 9.3346e-04 0.0011 0.0012 0.0013 0.0015 0.0016 0.0017 0.0019 0.0020 0.0021 0.0023 0.0024 0.0025 0.0027 0.0028 0.0029 0.0031 0.0032 0.0033 0.0035 0.0036 0.0037 0.0039 … ] y: [4×7500 double] yp: [4×7500 double] stats: [1×1 struct]
figure
yyaxis left
plot(sol.x, sol.y(1:3,:))
yyaxis right
plot(sol.x, sol.y(4,:))
grid
legend(compose('y_%d',1:size(sol.y,1)), 'Location','S')
function F = myfun(t, y, Pem,Peh, Da, b, U, yw, yf)
% global Pem Peh Da b U yw
y1 = y(1);
y2 = y(2);
y3 = y(3);
y4 = y(4);
F(1) = y2;
F(2) = Peh*(y2 + b*Da*y3*exp(-1/y1) - U*(yw-y1));
F(3) = y4;
F(4) = Pem*(y4 + Da*y3*exp(-1/y1));
F = F(:);
end
function f = mybc(ya,yb, Pem, Peh, yf)
% global Pem Peh yf
f(1) = ya(2)+Peh*(yf-ya(1));
f(2) = yb(2);
f(3) = ya(4) + Pem*(1-ya(3));
f(4) = yb(4);
f=f(:);
end
function g = guess(z) % initial guess for y and y'
g = [0.5
0.5
0.5
0.5];
end
Experiment to get different results.
.

More Answers (0)

Categories

Find more on Customize Object Indexing in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!