Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 3.23252e-20.

Hey!

I was trying to solve these equations and was getting an error as : Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 3.23252e-20. > In bvp4c (line 288) in compressiblechannel (line 15). How should I resolve it. here's my code and equations.

Thanks for your help in advance.

code:

close all
clear all
G = -1;
%Na = 0.86;
gamma = 1.4;
Pr = 0.71;
M = 0.1;
xstart = -1;
xend = 1;
nx = 51;
x = linspace(xstart,xend,nx);
solinit = bvpinit(x, [0;0;1;0]);
sol = bvp4c(@(y,z)bvpfcn(y,z,G, Pr, M, gamma), @bcfcn, solinit);

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 3.23252e-20.

% Extract variables from solution
U     = sol.y(1,:);  % U from solution
dUdy  = sol.y(2,:);  % dU/dy from solution
T0    = sol.y(3,:);  % T0 from solution
dT0dy = sol.y(4,:);  % dT0/dy from solution
y     = sol.x;       % y values (independent variable)
%Calculate d^2U/dy^2 using the equation from bvpfcn
% d2Udy2 = (dUdy .* dT0dy) + G * exp(T0);
% Plot d^2U/dy^2 vs y
figure(3)
% plot(d2Udy2, y);
% hold on
plot(diff(dUdy), y(2:end));
xlabel('d^2U/dy^2');
ylabel('y');
title('Plot of d^2U/dy^2 vs y');

figure(1)
plot(sol.y(1,:), sol.x)
xlabel ('U')
ylabel ('y')

figure(2)
plot(sol.y(3,:), sol.x)
xlabel ('T0')
ylabel ('y')

function dzdy = bvpfcn(y,z,G,Pr,M,gamma )
    U     = z(1);
    dUdy  = z(2);
    T0    = z(3);
    dT0dy = z(4);
    dzdy    = zeros(4,1);
    dzdy(1) = dUdy;
    dzdy(2) = -dUdy*dT0dy*0.76*(T0^ -1) + G/(T0^0.76);
    dzdy(3) = dT0dy;
    dzdy(4) = -Pr*M^2*(gamma-1)*((T0^0.76)*(dUdy^2));
end
function res = bcfcn(za,zb)
    res = [za(1);zb(1);za(3);zb(3)-1];
end

equation:

where mu0 = T0^ 0.76 and boundary conditions are: U = 0 at y= +1 and y = -1 and T0 = 0 at y =-1 and T0 = 1 at y = +1

[Torsten]

It seems mu0 = exp(-T0) was a better choice than mu0 = T0^0.76 for "bvp4c".

Since you divide by T0^0.76, T0 mustn't become zero anywhere. But you set T0 = 0 at y = -1. This might be a problem although it's only in a boundary point.

Another problem is that T0 should not become negative in the course of the iterations because T0^0.76 would produce a complex number.