How can I resolve the "Unable to solve collocation" error in MATLAB when trying to solve a nonlinear differential equation with nonlinear boundary conditions?

MIT6()

Pi = 1×251

0 0.010000000000000 0.020000000000000 0.030000000000000 0.040000000000000 0.050000000000000 0.060000000000000 0.070000000000000 0.080000000000000 0.090000000000000 0.100000000000000 0.110000000000000 0.120000000000000 0.130000000000000 0.140000000000000 0.150000000000000 0.160000000000000 0.170000000000000 0.180000000000000 0.190000000000000 0.200000000000000 0.210000000000000 0.220000000000000 0.230000000000000 0.240000000000000 0.250000000000000 0.260000000000000 0.270000000000000 0.280000000000000 0.290000000000000

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

Error in solution>MIT6 (line 26)
sol=bvp4c(@BVP_ODE,@BC,solinit)%options );

%% 
function  MIT6
clc
close all
format long
global A B eta
A =10 ;  % INTERNAL RADIUS 
B=200;   %% EXTERNEL RADIUS 
eta = B / A;  % THICKNESS RATIO
Mesh=100  ;       %% MESHING 



Pi = 0:0.01:2.5          %%  APPLIED INNER PRESSURE 

R = linspace(A, B, Mesh);
Lambda_theta = zeros(length(Pi), 1);
r_at_A = zeros(length(Pi), length(R));  %%% Specifying the size of matrix containing internal radius


for ik = 1:length(Pi)
    
    %options = bvpset('stats', 'on', 'NMax', 50000);
    solinit=bvpinit(R ,@Guess);
    sol=bvp4c(@BVP_ODE,@BC,solinit)%options );
    zz = deval(sol,R );
    r_at_A(ik,:) = zz(1,:);
    Lambda_theta(ik) = zz(1, 1) / A;
    display(r_at_A)
    figure(1)
    hold on
    plot(sol.x,sol.y(1,:),'LineWidth',3)
    
end

%%%%%%  deformed inner radius at A
a_at_A = (r_at_A(:,1))/A;





%% GRAPHICAL PART  %%%%% 
figure(2)
hold on
set(gca,'FontSize',16)
plot(a_at_A,Pi,'lineWidth',3)
set(get(gca,'Ylabel'),'Rotation',0)
ylabel('$P$','Interpreter','LaTeX','FontSize',20, 'FontWeight', 'normal', 'FontName', 'Times');
%         xlabel('$\Lambda_a ','Interpreter','LaTeX','FontSize',20, 'FontWeight', 'normal', 'FontName', 'Times');
xlabel('$\lambda_a$', 'Interpreter', 'latex', 'FontSize', 20, 'FontWeight', 'normal', 'FontName', 'Times');

grid on
hLegend=legend(['\eta =', num2str(eta)],'Location','SE')
% legend('boxoff')
% legend('Orientation','vertical')
set(hLegend, 'FontSize',18, 'Position', [0.7, 0.8, 0.1, 0.1]);
% Définir une boîte autour du graphe
ax = gca; % Récupérer l'objet axes actuel
set(ax, 'Box', 'on'); % Activer la boîte autour du graphe
%         title('$\alpha = 5$', 'FontSize', 16, 'FontWeight', 'bold', 'Color', 'black', 'Interpreter', 'latex');

%         %%
%          % Add annotations
%          [~, maxIndex] = max(Pi);
%          text(a_at_A(maxIndex), Pi(maxIndex), 'Maximum', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
%          xlabel('a/A');
%          ylabel('Pi');
%          title('Pi vs. a/A');
%         
%         
%%        
%         Ajouter la boîte de texte
annotation('textbox', [0.2, 0.6, 0.1, 0.1], 'String', ['Pmax : ', num2str(Pmax)],...
    'FitBoxToText', 'on', 'BackgroundColor','white', 'EdgeColor', 'black','Interpreter', 'latex');




%%
%         % save('pqfile.txt', 'a_at_A', '-ascii');

% Save Lambda_theta to a file (Excel format)
excelFileName = 'Lambda_theta.xlsx';                                
Lambda_theta_cell = num2cell(Lambda_theta);                                        
%xlswrite(excelFileName, Lambda_theta_cell);





%%    % STTING  THE SYSYETME OF ODE
    function dxdy = BVP_ODE(R, y)
        dxdy = [y(2);
            (2/3)*y(2)/R-(2/3)*(y(2)^4/y(1)^3)*R^2];
    end


%%
% SETING THE BOUNDARY CONDITION
    function res = BC(ya, yb)
        
        res = [1 - A^2 / (ya(1)^2 * ya(2)^3) + Pi(ik);
            1 - B^2 / (yb(1)^2 * yb(2)^3)];
    end


    function U = Guess(R)
        U = [  1  % Radial displacement profile
            1];                                       % Initial guess for the second variable
        
        % Initial guess for the second variable
    end
end 

[Star Strider]

The differential equations both increase exponentially and eventually go beyond 1E+40. Then it throws that error.

[Torsten]

Your boundary conditions don't fix y and y' at the endpoints of the R-interval.

Remember the infinite number of combinations for y and y' that give a certain value for y^2 * y'^3 !

[MOSLI KARIM]

sir Torsten can you offer me a solution please

[Torsten]

I don't know the background of your problem. Thus it's impossible to offer a "solution".

[MOSLI KARIM]

MIT6()

%% 
function  MIT6
    clc
    close all
    format long
    global A B eta
    A =10 ;  % INTERNAL RADIUS 
    B=200;   %% EXTERNEL RADIUS 
    eta = B / A;  % THICKNESS RATIO
    Mesh=100  ;       %% MESHING 
    Pmin =0;
    Pmax=2.5;
   
    N =100;
     %h = (Pmax - Pmin) /N;
     Pi = 0:2.5;       %%  APPLIED INNER PRESSURE 
%     Pi=linspace(0,2.51,20  )
    R = linspace(A, B, Mesh);
    Lambda_theta = zeros(length(Pi), 1);
    r_at_A = zeros(length(Pi), length(R));  %%% Specifying the size of matrix containing internal radius
  
    
    for ik = 1:length(Pi)
 
        %options = bvpset('stats', 'on', 'NMax', 50000);
        solinit=bvpinit(R ,@Guess);
        sol=bvp4c(@BVP_ODE,@BC,solinit);%options );
        zz = deval(sol,R );
        r_at_A(ik,:) = zz(1,:);
        Lambda_theta(ik) = zz(1, 1) / A;
        %display(r_at_A)
        figure(1)
        hold on
        plot(sol.x,sol.y(1,:),'LineWidth',3)
        figure(2)
        Ur=sol.y(1,:)-sol.x;
        figure(2)
        hold on 
        plot(sol.x,Ur)
    end

        %%%%%%  deformed inner radius at A
        a_at_A = (r_at_A(:,1))/A;
   



        
        %% GRAPHICAL PART  %%%%% 
        figure(3)
        hold on
        set(gca,'FontSize',16)
        plot(a_at_A,Pi,'lineWidth',3)
        set(get(gca,'Ylabel'),'Rotation',0)
        ylabel('$P$','Interpreter','LaTeX','FontSize',20, 'FontWeight', 'normal', 'FontName', 'Times');
%         xlabel('$\Lambda_a ','Interpreter','LaTeX','FontSize',20, 'FontWeight', 'normal', 'FontName', 'Times');
        xlabel('$\lambda_a$', 'Interpreter', 'latex', 'FontSize', 20, 'FontWeight', 'normal', 'FontName', 'Times');

        grid on
        hLegend=legend(['\eta =', num2str(eta)],'Location','SE');
        % legend('boxoff')
        % legend('Orientation','vertical')
        set(hLegend, 'FontSize',18, 'Position', [0.7, 0.8, 0.1, 0.1]);
        % Définir une boîte autour du graphe
        ax = gca; % Récupérer l'objet axes actuel
        set(ax, 'Box', 'on'); % Activer la boîte autour du graphe
%         title('$\alpha = 5$', 'FontSize', 16, 'FontWeight', 'bold', 'Color', 'black', 'Interpreter', 'latex');

%         %%
%          % Add annotations
%          [~, maxIndex] = max(Pi);
%          text(a_at_A(maxIndex), Pi(maxIndex), 'Maximum', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
%          xlabel('a/A');
%          ylabel('Pi');
%          title('Pi vs. a/A');
%         
%         
%%        
%         Ajouter la boîte de texte
annotation('textbox', [0.2, 0.6, 0.1, 0.1], 'String', ['Pmax : ', num2str(Pmax)],...
           'FitBoxToText', 'on', 'BackgroundColor','white', 'EdgeColor', 'black','Interpreter', 'latex');
        
      
        
        
        %%
%         % save('pqfile.txt', 'a_at_A', '-ascii');
        
        % Save Lambda_theta to a file (Excel format)
        %excelFileName = 'Lambda_theta.xlsx';                                
        %Lambda_theta_cell = num2cell(Lambda_theta);                                        
        %xlswrite(excelFileName, Lambda_theta_cell);
 


        
        
        %%    % STTING  THE SYSYETME OF ODE
        function dxdy = BVP_ODE(R, y)
            dxdy = [y(2);
                 (2/3)*y(2)/R-(2/3)*(y(2)^4/y(1)^3)*R^2];
        end

     
%%
        % SETING THE BOUNDARY CONDITION
        function res = BC(ya, yb)

            res = [1 - A^2 / (ya(1)^2 * ya(2)^3) + Pi(ik);
                1 - B^2 / (yb(1)^2 * yb(2)^3)];
        end
 
   function U = Guess(R)
     U = [ R  ; % Radial displacement profile
        1];                                       % Initial guess for the second variable
 
                          
   end

    
end

 

This MATLAB code describes the inflation phenomenon of a hollow sphere with an inner radius A and an outer

radius B, composed of a compressible hyperelastic material of Blatz-Ko type. My issue lies in the fact that the pressure curve as a function of the deformed inner radius does not decrease once it reaches its maximum value."

[MOSLI KARIM]

this my equation sir

[Torsten]

My issue lies in the fact that the pressure curve as a function of the deformed inner radius does not decrease once it reaches its maximum value.

Looking at your solution curves, I don't know what you mean (see above).