How to integrate a function over whose integrand upper limit is defined as an array?

Suppose I have this equation.

I am trying to graph η vs. θ for . I have a table of at given η values such as

I tried integration using arrays but that proved difficult. So then I plotted the graph of and came out with best fit curves, so my thiniking was to simply plug the function into the above equation. But that also proved difficult. How can I accomplish the integration of the above eqaution and then plot the results? Any guidance helps.

I only have this

syms eta Pr
eta = 0:0.1:8;
f_d = -1*10.^-5*eta.^6 + 0.0002*eta.^5 + 0.0012*eta.^4 - 0.0196*eta.^3 + 0.0345*eta.^2 + 0.3118*eta + 0.002;
f_dd = 6*10^-5*eta.^6 - 0.0015*eta.^5 + 0.0129*eta.^4 - 0.0458*eta.^3 + 0.0428*eta.^2 - 0.0177*eta + 0.3332;

eta_array = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.6 3 3.4 3.8 4.2 4.6 5 5.6 6.2 7 8];
f_dd_array = [0.33206 0.33199 0.33147 0.33008 0.32739 0.32301 0.31659 0.30787 0.29667 0.28293 0.26675 0.24835 0.20646 0.16136 0.11788 0.08013 0.05052 0.02948 0.01591 0.00543 0.00155 0.00022 0.00001];

intEta = int(f_dd,eta,[0 eta]); %numerator
intInf = int(f_dd,eta,[0 inf]); %denomenator

theta = 1 - intEta/intInf;

Answer 1

[Alan Stevens]

You could try the following

eta_array = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.6 3 3.4 3.8 4.2 ...
             4.6 5 5.6 6.2 7 8];
f_dd_array = [0.33206 0.33199 0.33147 0.33008 0.32739 0.32301 0.31659 0.30787 ...
    0.29667 0.28293 0.26675 0.24835 0.20646 0.16136 0.11788 0.08013 0.05052 ...
    0.02948 0.01591 0.00543 0.00155 0.00022 0.00001];

Pr = 30;

den = trapz(eta_array,f_dd_array.^Pr);  % denominator 
disp(den)
num = zeros(1,numel(eta_array));

for i = 2:numel(eta_array)
    
    num(i) = trapz(eta_array(1:i),f_dd_array(1:i).^Pr);
     
end
 theta = 1 - num/den;
 
 plot(eta_array,theta,'-o'),grid
 xlabel('\eta'), ylabel(['\theta for Pr = ',num2str(Pr)])

This works well enough for values of Pr up to 30, but raising the f'' values to the power of 1000 results in an array of what are essentially zeros!

[Missael Hernandez]

Wow this is very clear. Thank you! How would I loop through multiple Pr values?

[Alan Stevens]

Something like this:

eta_array = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.6 3 3.4 3.8 4.2 ...
             4.6 5 5.6 6.2 7 8];
f_dd_array = [0.33206 0.33199 0.33147 0.33008 0.32739 0.32301 0.31659 0.30787 ...
    0.29667 0.28293 0.26675 0.24835 0.20646 0.16136 0.11788 0.08013 0.05052 ...
    0.02948 0.01591 0.00543 0.00155 0.00022 0.00001];

Pr = [0.5, 1, 10, 20, 30];

for j = 1:numel(Pr)
    den = trapz(eta_array,f_dd_array.^Pr(j));  % denominator 
    %disp(den)
    num = zeros(1,numel(eta_array));

    for i = 2:numel(eta_array)

        num(i) = trapz(eta_array(1:i),f_dd_array(1:i).^Pr(j));

    end
     theta(j,:) = 1 - num/den;

end
 
 plot(eta_array,theta,'-o'),grid
 xlabel('\eta'), ylabel(['\theta '])

Answer 2

[Walter Roberson]

Your posted is a polynomial that is negative between approximately 4 and 11.

When Pr is not an integer, such as then the denominator integral could be considered to be infinite plus a finite complex component. The complex component of the numerator integral depends upon the end point of the integration, which is variable in your system.

However... No matter whether the numerator is pure finite real or is finite real plus finite complex, by multiplying and dividing by conj() of the denominator, and expanding, and doing some quick infinite limits, you can show that the ratio is going to come out as 0 (not real 0 and some complex component.)

So you have theta being 1 minus 0, so theta is going to be 1.

There are potential exceptions at boundary points such as or

[Alan Stevens]

But the polynomial is clearly a poor fit to the data

[Walter Roberson]

format long g
eta_array = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.6 3 3.4 3.8 4.2 ...
             4.6 5 5.6 6.2 7 8];
f_dd_array = [0.33206 0.33199 0.33147 0.33008 0.32739 0.32301 0.31659 0.30787 ...
    0.29667 0.28293 0.26675 0.24835 0.20646 0.16136 0.11788 0.08013 0.05052 ...
    0.02948 0.01591 0.00543 0.00155 0.00022 0.00001];
syms eta
f_dd = 6*10^-5*eta.^6 - 0.0015*eta.^5 + 0.0129*eta.^4 - 0.0458*eta.^3 + 0.0428*eta.^2 - 0.0177*eta + 0.3332;
p = polyfit(eta_array, f_dd_array, 6);
mat2str(p)

ans = '[5.87443024024372e-05 -0.00145347836835335 0.01287615621175 -0.0457872748212294 0.0428018872107394 -0.0177415790268976 0.333171626268646]'

vpa(p,3)

ans = 

sym_f_dd = poly2sym(p, eta);
scatter(eta_array, f_dd_array, 'bo')
hold on
fplot(f_dd, [0 8], 'k-');
fplot(sym_f_dd, [0 8], 'g-');
legend({'data', 'f_dd given', 'polyfit'})

You can see that the coefficients posted are more or less 3 digit round-offs of the polyfit() values. The round-off is enough to give a poor approximation at higher η but if you use higher precision then the fit is not bad

rp = roots(p);
rp(imag(rp) == 0)

ans = 2×1

7.99655276703378 6.93633138446102

●

And that shows that between about 6.9 and 7.99 the polynomial is negative. In such a case, when is not integral, complex integrals are going to result when taken over any finite upper bound that is greater than the 6.9-whatever lower bound shown here.

On the other hand, as I already explored... it does not matter to the end result if there is an imaginary component to the upper or lower integrals: the real component of the lower integral will be infinity, and that means the ratio of the two integrals will be 0.

[Missael Hernandez]

the fitted curve i gave, was an equation Excel derived. that is weird that was the outcome. It seem to fit the scatter points well enough. Perhaps I didn't coupy-and-paste correctly. thank you for the input. it is greatly appreciated.