![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/832125/image.png)
Coupled Mass and energy equations using bvp4c
4 views (last 30 days)
Show older comments
Jagadeesh Korukonda
on 12 Dec 2021
Commented: Star Strider
on 12 Dec 2021
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/832085/image.png)
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);
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
0 Comments
Accepted Answer
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.)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/832125/image.png)
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)
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.
.
2 Comments
More Answers (0)
See Also
Categories
Find more on Customize Object Indexing 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!