heat and mass transfer problem nonlinear differential equations

The code is running, the graph is also came. But it does not give the convergence for the curves s, , θ, and f.

fluid_Sol_bvp4c2

Warning: Unable to meet the tolerance without using more than 1428 mesh points.
The last mesh of 883 points and the solution are available in the output argument.
The maximum residual is 75.7887, while requested accuracy is 0.001.

y = 7×883

0 -0.0007 -0.0029 -0.0065 -0.0115 -0.0157 -0.0205 -0.0260 -0.0321 -0.0365 -0.0412 -0.0462 -0.0515 -0.0570 -0.0628 -0.0690 -0.0753 -0.0820 -0.0890 -0.0962 -0.1037 -0.1115 -0.1196 -0.1279 -0.1365 -0.1454 -0.1545 -0.1640 -0.1737 -0.1837 0 -0.0570 -0.1142 -0.1715 -0.2288 -0.2669 -0.3051 -0.3431 -0.3811 -0.4064 -0.4317 -0.4569 -0.4820 -0.5071 -0.5322 -0.5572 -0.5821 -0.6070 -0.6318 -0.6566 -0.6813 -0.7060 -0.7305 -0.7551 -0.7795 -0.8039 -0.8282 -0.8525 -0.8766 -0.9007 -2.2514 -2.2615 -2.2672 -2.2693 -2.2683 -2.2662 -2.2632 -2.2593 -2.2546 -2.2511 -2.2474 -2.2433 -2.2391 -2.2346 -2.2300 -2.2251 -2.2201 -2.2150 -2.2097 -2.2043 -2.1987 -2.1931 -2.1873 -2.1815 -2.1756 -2.1696 -2.1635 -2.1573 -2.1511 -2.1449 1.0000 1.0021 1.0039 1.0055 1.0069 1.0078 1.0086 1.0094 1.0102 1.0107 1.0112 1.0116 1.0121 1.0126 1.0130 1.0135 1.0139 1.0143 1.0148 1.0152 1.0156 1.0160 1.0164 1.0169 1.0173 1.0177 1.0181 1.0185 1.0189 1.0193 0.0894 0.0761 0.0665 0.0594 0.0540 0.0511 0.0486 0.0466 0.0448 0.0438 0.0429 0.0420 0.0413 0.0406 0.0400 0.0394 0.0389 0.0384 0.0380 0.0376 0.0373 0.0370 0.0367 0.0364 0.0362 0.0360 0.0358 0.0356 0.0354 0.0353 1.0000 1.3878 1.7090 1.9750 2.1952 2.3205 2.4310 2.5285 2.6144 2.6659 2.7133 2.7569 2.7970 2.8338 2.8677 2.8989 2.9275 2.9538 2.9780 3.0003 3.0208 3.0396 3.0569 3.0728 3.0874 3.1009 3.1132 3.1246 3.1351 3.1447 16.8518 13.9562 11.5569 9.5694 7.9232 6.9861 6.1597 5.4310 4.7884 4.4029 4.0483 3.7223 3.4225 3.1469 2.8934 2.6604 2.4461 2.2490 2.0679 1.9013 1.7482 1.6073 1.4778 1.3588 1.2493 1.1487 1.0561 0.9710 0.8928 0.8209

function fluid_Sol_bvp4c2
    Pr = 10;
    Le = 10;
    Nr = 0.5;
    Nb = 0.5;
    Nt = 0.5;
    
    % Defining parameters 
    solinit = bvpinit(linspace(0, 10, 100), [0 0 1 1 0 1 0]);
    sol     = bvp4c(@bvp2D, @bc2D, solinit);
    
    x = sol.x;
    y = sol.y
    
    % Plotting of the functions
    figure(1)
    plot(x, y(1,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('s(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of the velocity
    figure(2)
    plot(x, y(2,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('s^{\prime}(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of theta
    figure(3)
    plot(x, y(4,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta ', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('\theta ', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of the fun
    figure(4)
    plot(x, y(6,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('f(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off

% Residual of the boundary conditions
    function res = bc2D(y0, yinf)
        res = [y0(1);           % y1 to start at 0
               y0(2);           % y2 to start at 0
               y0(4) - 1;       % y4 to start at 1
               y0(6) - 1;       % y6 to start at 1
               yinf(2) - 0;     % desire y2 to converge to 0
               yinf(4) - 0;     % desire y4 to converge to 0
               yinf(6) - 0];    % desire y5 to converge to 0
    end

% System of First Order ODEs
    function dydx = bvp2D(~, y)
        yy1  = - 1/(4*Pr)*(3*y(1)*y(3) - 2*y(2)^2) - y(4) + Nr*y(6);
        yy2  = - 3/4*y(1)*y(5) - Nb*y(7)*y(5) - Nt*y(5)^2;
        yy3  = - 3/4*Le*y(7)   - (Nt/Nb)*yy2;
        dydx = [y(2); y(3); yy1; y(5); yy2; y(7); yy3];
    end
end

I dont know whether i have to use extend code in bvp4c. Please help me do the needful. Thanks in advance.

[Torsten]

Please include equations and boundary conditions in a mathematical notation so that we can check whether your implementation is correct.

Answer 1

[Sam Chak]

Hi @uma

I wanted to inform you that I have discovered five sets of Equilibrium Points for the system. However, it appears that the system is unstable. To address this instability, it would be advisable to introduce some manipulatable variables into the system. As a temporary solution, I have incorporated a guess() function into your code. This adjustment has resulted in a different solution when using bvp4c().

%% Equilibrium Points
ye1 = zeros(1, 7);            
ye2 = [1, 0, 0, 0,   0, 0,   0];
ye3 = [0, 0, 0, 0.2, 0, 0.4, 0];
ye4 = [0, 0, 0, 0.4, 0, 0.8, 0];
ye5 = [0, 0, 0, 0.6, 0, 1.2, 0];

fluid_Sol_bvp4c2

Warning: Unable to meet the tolerance without using more than 1428 mesh points.
The last mesh of 714 points and the solution are available in the output argument.
The maximum residual is 1.6282, while requested accuracy is 0.001.

function fluid_Sol_bvp4c2
    Pr = 10;
    Le = 10;
    Nr = 0.5;
    Nb = 0.5;
    Nt = 0.5;
    
    % Defining parameters 
    a       = 0;
    b       = 10;
    xmesh   = linspace(a, b, 101);
    solinit = bvpinit(xmesh, @guess);
    sol     = bvp4c(@bvp2D, @bc2D, solinit);
    
    % solinit = bvpinit(linspace(0, 10, 100), [0 0 1 1 0 1 0]);
    % sol     = bvp4c(@bvp2D, @bc2D, solinit);
    
    x = sol.x;
    y = sol.y;
    
    % Plotting of the functions
    figure(1)
    plot(x, y(1,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('s(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of the velocity
    figure(2)
    plot(x, y(2,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('s^{\prime}(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of theta
    figure(3)
    plot(x, y(4,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta ', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('\theta ', 'fontweight', 'bold', 'fontsize', 16)
    hold off
    
    % Plotting of the fun
    figure(4)
    plot(x, y(6,:), 'linewidth', 1), grid
    hold on
    xlabel('\eta', 'fontweight', 'bold', 'fontsize', 16)
    ylabel('f(\eta)', 'fontweight', 'bold', 'fontsize', 16)
    hold off

% Residual of the boundary conditions
    function res = bc2D(y0, yinf)
        res = [y0(1);           % y1 to start at 0
               y0(2);           % y2 to start at 0
               y0(4) - 1;       % y4 to start at 1
               y0(6) - 1;       % y6 to start at 1
               yinf(2) - 0;     % desire y2 to converge to 0
               yinf(4) - 0;     % desire y4 to converge to 0
               yinf(6) - 0];    % desire y5 to converge to 0
    end

% System of First Order ODEs
    function dydx = bvp2D(~, y)
        yy1  = - 1/(4*Pr)*(3*y(1)*y(3) - 2*y(2)^2) - y(4) + Nr*y(6);
        yy2  = - 3/4*y(1)*y(5) - Nb*y(7)*y(5) - Nt*y(5)^2;
        yy3  = - 3/4*Le*y(7)   - (Nt/Nb)*yy2;
        dydx = [y(2); y(3); yy1; y(5); yy2; y(7); yy3];
    end
end

% initial guess
function g = guess(x)
    g     = [sin(x);
             cos(x);
            -sin(x);
            -cos(x);
             sin(x);
             cos(x);
            -sin(x)];
end

[uma]

No Sam sir, I have to get the convergence of s from zero to upwards.. S' from zero to 10.theta from 1 to 4..f from 1 to 2 almost...

Answer 2

[Fareeha]

@uma your code seems to be okay. you just have to adjust parameters values and you will get fine results. Also better to use

yinf(2); yinf(4); yinf(6)]; instead yinf(2) - 0; yinf(4) - 0; yinf(6) - 0];. i hope this will help.