How do I input a lot of values and record each single output in Bvp4c?

Hello everyone,

In my code below, how can I make the code input the values of "Pe" from { 25 to 1025 } and record the generated output values for EACH single input number.

For example; when I input a value of "Pe" of 25, the generated value of "MainpADM" is { 0.0054 }, and another input value of "Pe" of { 26 } will generate { 0.0065 }. and of course I am recording each corresponding output

This old-fashioned way will take A LOT of time to reach the final 1025 input. Therefore, Could you kindly help me in this regard?

The three codes of bvp4c:

MainADM

global Pe Taw Yx B A; 
% 
 Pe = 61;
 Taw=4.5; % dimensionless residemnce time 
 Yx=0.212;
 B=0.482;
 A=1.01; % A=1+(Xo/So*Yx)
 xlow = 0;
 xhigh = 1;
 y=[];
 solinit = bvpinit(linspace(xlow,xhigh,20),[0 1]);
 sol = bvp4c(@bvpode5,@bvpbc5,solinit);
 plot(sol.x,sol.y(1,:),'r-')
 xlabel('Rxtr Length Z/L')
 ylabel('Concentration Fraction (S/So)')
 axis ([0 1 0 1]);
 grid on
%
legend('alpha')
fh=figure(1);
set(fh,'color','white')
title('Concentration Profiles Across the Rxtr')
disp (sol.y(1,end)); 

bvpode5

function dydx = bvpode5(x,y)
%
global Pe Taw Yx B A; 
%
dydx = [y(2) 
    Pe*y(2)+((Pe*Taw*y(1)*(A-y(1)))/(B*Yx*(A-y(1))+y(1)))];
end

bvpbc5

function res = bvpbc5(ya,yb)
global Pe Taw Yx B; 
%
res = [ ya(1)-1-((1/Pe)*ya(2))
    yb(2)];
end 

Answer 1

[Tejas]

Hi Hussein,

To record the values of 'MainpADM' for different values of 'Pe', follow these steps:

• Preallocate the variable 'MainpADM_results' to store the results of 'MainpADM' for each 'Pe' value.

    % MainADM
    global Pe Taw Yx B A; 
    Taw = 4.5;  
    Yx = 0.212;
    B = 0.482;
    A = 1.01; 
    % Range of Pe values
    Pe_values = 25:1025;
    % Initialize storage for results
    MainpADM_results = zeros(size(Pe_values));

• Use a 'for' loop to iterate over the values of 'Pe'. More information on 'for' loops can be found in this documentation: https://www.mathworks.com/help/matlab/ref/for.html .

    % Loop over all Pe values
    for i = 1:length(Pe_values)
        Pe = Pe_values(i);
        
        xlow = 0;
        xhigh = 1;
        solinit = bvpinit(linspace(xlow, xhigh, 20), [0 1]);
        sol = bvp4c(@bvpode5, @bvpbc5, solinit);
        
        % Record the result for the current Pe
        MainpADM_results(i) = sol.y(1, end);
    end

Answer 2

[Walter Roberson]

This old-fashioned way will take A LOT of time to reach the final 1025 input.

Looping is the main practical way to solve this problem. It can take a long time, and that is just how it is.

The less practical way to solve the problem is to transform it into a system of 2 * (1025 - 25 + 1) = 2 * 1001 variables, and vectorize your bvpode5 and bvpbc5. For example,

function dydx = bvpode5(x,y)
%
global Pe Taw Yx B A; 
%
dydx = [y(2:2:end), 
       Pe.*y(2:2:end) + ((Pe.*Taw.*y(1:2:end).*(A-y(1:2:end)))./(B.*Yx.*(A-y(1:2:end))+y(1:2:end)))];
dydx = dydx(:);
end

and then make a single call to bvp4c with Pe set to 25:1025. This will return all the solutions simultaneously