Coupled Mass and energy equations using bvp4c

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

Answer 1

[Star Strider]

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.

.

[Jagadeesh Korukonda]

Thank you so much

[Star Strider]

As always, my pleasure!

.