Solution of a system of ODEs by means of bvp4c

Hi eveyone,

I am trying to solve the following problem thorugh the bvp4c solver:

clear,clc,close all

%% coordinate
R=1e-3;

%% xmesh command
xmesh = linspace(1e-8,R,50);

%% bvp4c 
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@bvpfcn, @bcfcn, solinit);
plot(sol.x,sol.y(1,:))

where, in separated .m files, I defined:

function dydx = bvpfcn(r,y)

k_0=3508;
n=0.2;
T0=373.15;
a=0.01;
k=0.18;
k0=3508;
dP=100;
dZ=1e-3;

dydx = zeros(5,1);
dydx = [y(2)
    (k0*exp(a*T0-y(1))*(((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((n+1)/n))*(r/k)-y(2))*(1/r)
    (((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((1)/n))
    y(5)
    -(a*k0*exp(a*T0-y(1))*y(4)*(((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((n+1)/n))*(r/k)-y(5))*(1/r)];
end



function res = bcfcn(ya,yb)
Tw=500;
res = [ yb(3)
    yb(1)-Tw
    ya(2)
    yb(4)
    ya(5)];
end



function g = guess(x)
Tw=500;
g = [Tw
    x^2
    x^2
    x
    x^0.5];
end

The following error is returned:

Error using bvp4c (line 197)

Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in untitled5 (line 11)

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

I would appreciate Your help in finding the solution to this problem.

Best regards,

AP

Answer 1

[Torsten]

Works if you start at r=1e-8 instead of r=0 (see above).

Answer 2

[Star Strider]

The singlular Jacobian was likely the result of the ‘(1/r)’ term, since ‘r’ initially begins at 0. While changing that removed the error, it likely does not do much for the result —

%% coordinate
R=1e-3;
%% xmesh command
xmesh = linspace(1E-8,R,50);
%% bvp4c 
solinit = bvpinit(xmesh, @guess);
sol = bvp4c(@bvpfcn, @bcfcn, solinit)

sol = struct with fields:
    solver: 'bvp4c'
         x: [1.0000e-08 2.0418e-05 4.0826e-05 6.1234e-05 8.1642e-05 1.0205e-04 1.2246e-04 1.4287e-04 1.6327e-04 1.8368e-04 2.0409e-04 2.2450e-04 2.4491e-04 2.6531e-04 2.8572e-04 3.0613e-04 3.2654e-04 3.4695e-04 3.6735e-04 3.8776e-04 4.0817e-04 … ]
         y: [5×50 double]
        yp: [5×50 double]
     stats: [1×1 struct]

r = sol.x;
y = sol.y;

figure
semilogy(r, real(y))
grid

figure
semilogy(r, abs(y))
grid

Warning: Negative data ignored

% where, in separated .m files, I defined:
function dydx = bvpfcn(r,y)
k_0=3508;
n=0.2;
T0=373.15;
a=0.01;
k=0.18;
k0=3508;
dP=100;
dZ=1e-3;
dydx = zeros(5,1);
dydx = [y(2)
        (k0*exp(a*T0-y(1))*(((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((n+1)/n))*(r/k)-y(2))*(1/r)
        (((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((1)/n))
        y(5)
        -(a*k0*exp(a*T0-y(1))*y(4)*(((dP/dZ)*(r/2)*(1/k0*exp(a*(T0-y(1)))))^((n+1)/n))*(r/k)-y(5))*(1/r)];
end
function res = bcfcn(ya,yb)
Tw=500;
res = [ yb(3)
        yb(1)-Tw
        ya(2)
        yb(4)
        ya(5)];
end
function g = guess(x)
Tw=500;
g = [Tw
     x^2
     x^2
     x
     x^0.5];
end

.