Using 'for' loop in integral correctly

Hello! I need to solve an integral which is definite and has an unknown parameter called 'B'. Its limits are t1(i) and t2(i). Because of i=1:3, I need to use the for loop. I'm trying to do it, but it's giving me results with long numbers and I felt like I'm doing something wrong. I'm sending the picture of integral I suppose to solve.

And here is the code in matlab I wrote for that integral:

for i=1:3

syms B t

y(i)=int((1-(B/t)+(1/2)*(B^2/t^2)-(1/6)*(B^3/t^3))*(a(i)+(b(i)*t)+(c(i)*t^2)),t1(i),t2(i))

end

y=y'

I tried to have the results as 'y(i)'. Because I will need yi(1) and yi(2) values for determining 'B' at the next step. I am expecting your help.

Answer 1

[Alan Stevens]

Just define y as a function of t, B, a, b, and c, wihout regard to i. Then call the function with the appropriate parameters when needed. ie. something like:

y = @(t,B,a,b,c) int(1- B/t +(1/2)*(B/t)^2 + ... etc.);

[Alan Stevens]

On the other hand it might be simpler just to do the integral once and then define the function, something like the following example:

yfn = @(t,a,b,c,B) t^2*(b/2 - (B*c)/2) - log(t)*((c*B^3)/6 - (b*B^2)/2 + a*B) + (c*t^3)/3 + ((B^3*a)/2 - t*(- b*B^3 + 3*a*B^2))/(6*t^2) + t*((c*B^2)/2 - b*B + a);

a =1; b = 2; c = 3;
t1 = 1;
t2 = 5;
B = 0.5;
deltay = yfn(t2,a,b,c,B) - yfn(t1,a,b,c,B)

[Enver Can Kılıç]

Thank you for your answer, but I still can't do it. Can you give me in more detailed way? I am so new in Matlab. By the way, the parameters a(i), b(i), c(i), t1(i) and t2(i) are known. I found them at previous steps. The only unknown parameter is B.

[Alan Stevens]

Here is one way to do it (somewhat clumsy, but it works!). I've assumed you are trying to find the values of B that make the integral match certain target values:

% Arbitrary data (replace with your own known data
a = [1 2 3];
b = [0.5 0.6 0.7];
c = [0.5 1.5 2.5];

t1 = [1 2 3];
t2 = [5 6 7];

ytarget = [2 3 4];


B0 = 1;  % Initial guess 

for i = 1:3  % Call fzero 3 times to find the three values of B
    
    B(i) = fzero(@yfn,B0,[],t1(i),t2(i),a(i),b(i),c(i),ytarget(i));
    
end

disp(B)

% The integral in the following was obtained with symbolic toolbox 
function yzero = yfn(B,t1,t2,a,b,c,ytarget)

         y2 = t2^2*(b/2 - (B*c)/2) - log(t2)*((c*B^3)/6 - (b*B^2)/2 + a*B) + (c*t2^3)/3 + ((B^3*a)/2 - t2*(- b*B^3 + 3*a*B^2))/(6*t2^2) + t2*((c*B^2)/2 - b*B + a);
         y1 = t1^2*(b/2 - (B*c)/2) - log(t1)*((c*B^3)/6 - (b*B^2)/2 + a*B) + (c*t1^3)/3 + ((B^3*a)/2 - t1*(- b*B^3 + 3*a*B^2))/(6*t1^2) + t1*((c*B^2)/2 - b*B + a);
         yzero = y2 - y1 - ytarget;
end

[Enver Can Kılıç]

Dear Alan,

The B value has to be found at the next step after the integral. And it is only one value, it should be found 0.06225. I did the integral once as you said and then tried to get three results of integral. I only need first two results of integral to find the B at the next step. I called that first two values y1 and y2. I did something like this after your answers:

syms t a b c B

y=int((1-(B/t)+(1/2)*(B^2/t^2)-(1/6)*(B^3/t^3))*(a+(b*t)+(c*t^2)))

yfn=@(t,a,b,c,B) t^2*(b/2 - (B*c)/2) - log(t)*((c*B^3)/6 - (b*B^2)/2 + a*B) + (c*t^3)/3 + ((B^3*a)/2 - t*(- b*B^3 + 3*a*B^2))/(6*t^2) + t*((c*B^2)/2 - b*B + a);

a=1.223570649023627; b =-11.041168221588103; c =23.521235780695218;

t1=0.021206896551724;

t2=0.135159216180759;

y1=yfn(t2,a,b,c,B)-yfn(t1,a,b,c,B)

yfn=@(t,a,b,c,B) t^2*(b/2 - (B*c)/2) - log(t)*((c*B^3)/6 - (b*B^2)/2 + a*B) + (c*t^3)/3 + ((B^3*a)/2 - t*(- b*B^3 + 3*a*B^2))/(6*t^2) + t*((c*B^2)/2 - b*B + a);

a=1.201007207988730; b =-9.767715978553012; c =19.356045435699443;

t1=0.021494252873563;

t2=0.177302497494984;

y2=yfn(t2,a,b,c,B)-yfn(t1,a,b,c,B)

syms B

denklem=y1/y2==k/h

(k and h are found at previous step, y1 and y2 are first two results of the integral with the unknown 'B'. y3 is not needed to find the 'B' value at this work.)

cevap=solve(denklem,B)

But it gave me three B results and the results are starting with 'root'. I am sending you the work that I try to write at matlab. The integral is at page 3. And at page 4, B value is found 0.06225.

Thank you so much for your attention and help. I am just stucked at this integral stage. If I can solve this, the rest of the work is easy to solve.

Best regards,

Can

[Alan Stevens]

Actually, the result in sumodelson2000 is B = 0.062973573741391. B is then redefined as 0.06225 for the subsequent calculations. It looks like that paper used MathCad to do the calculations. I've copied and pasted the values of their constants in the program below. However, not all were presented with high precision, so I suspect that is causing the Matlab result for B to be slightly different.

I've done the integration (as an indefinite integral) separately, and coded the result into Matlab - see below.

The following results in B = 0.0731.

B0 = 0.1;  % Initial guess
B = fzero(@rfn,B0);
disp(B)

function r = rfn(B)

        a = [1.2235706, 1.2010072, 1.188249]; 
        b = [-11.0411682, -9.767716, -9.0103308];
        c = [23.5212358, 19.3560454, 16.8841345];

        t1 = [0.0212069, 0.0214943, 0.0217816];
        t2 = [0.1351592, 0.1773025, 0.2120814];

        p = [1.6586, 3.6583, 4.514]*10^-5;

        rho = [0.9956, 0.9922, 0.9880]*10^6;
        T = [303, 313, 323];

        betaT = [0.2688823, 0.271018, 0.2736945];

        k = p(1)/p(2);
        h1 = rho(1)*T(1)/betaT(1);
        h2 = rho(2)*T(2)/betaT(2);
        h = h1/h2;
        DI1 = yfnintegral(t2(1),B,a(1),b(1),c(1)) - yfnintegral(t1(1),B,a(1),b(1),c(1));
        DI2 = yfnintegral(t2(2),B,a(2),b(2),c(2) ) - yfnintegral(t1(2),B,a(2),b(2),c(2) );
        r = h*DI1/DI2 - k;
         
end


function I = yfnintegral(t,B,a,b,c)

         I1 = t.^2*(b - B*c)/2;
         I2 = -log(t)*(a*B - b*B^2/2 + c*B^3/6);
         I3 = c*t.^3/3;
         I4 = -( t*(a*B^2/2 -b*B^3/6) - a*B^3/12 )./t.^2;
         I5 = t*(a - b*B + c*B^2/2);
         I = I1 + I2 + I3 + I4 + I5;
end

[Enver Can Kılıç]

That difference is not important, thanks! Did you do the integral in matlab or somewhere else? Because integral result is different from matlab's in here. Now I try to code the next step to find 'Adi' values. But it gives an error like this: 'Function definitions in a script must appear at the end of the file.

Move all statements after the "yfnintegral" function definition to before the

first local function definition.'

Should I try to code last two steps as functions for avoiding this error?

[Alan Stevens]

I did the integration in MathCad - it's easier to see the result - not (quite) as many brackets involved!

You could simply keep moving the yfnintegral function definition down until it appears after all your other calculations.

[Enver Can Kılıç]

When I move the yfnintegral function definiton down and try to do next other calculations before than it, it can't take into account some parameters like DI1 and DI2. At function steps some parameters (DI1 and DI2) are defined and I need them to be taken into account to do next steps.

[Alan Stevens]

Ok, here's a version that allows you to continue the calculations easily:

a = [1.2235706, 1.2010072, 1.188249]; 
b = [-11.0411682, -9.767716, -9.0103308];
c = [23.5212358, 19.3560454, 16.8841345];

t1 = [0.0212069, 0.0214943, 0.0217816];
t2 = [0.1351592, 0.1773025, 0.2120814];

p = [1.6586, 3.6583, 4.514]*10^-5;

rho = [0.9956, 0.9922, 0.9880]*10^6;
T = [303, 313, 323];

betaT = [0.2688823, 0.271018, 0.2736945];

k = p(1)/p(2);
h1 = rho(1)*T(1)/betaT(1);
h2 = rho(2)*T(2)/betaT(2);
h = h1/h2;

B0 = 0.1;  % Initial guess
B = fzero(@rfn,B0,[],a,b,c,t1,t2,k,h);
disp(B)

DI1 = yfnintegral(t2(1),B,a(1),b(1),c(1)) - yfnintegral(t1(1),B,a(1),b(1),c(1));
DI2 = yfnintegral(t2(2),B,a(2),b(2),c(2) ) - yfnintegral(t1(2),B,a(2),b(2),c(2) );

% Insert the other calculations here before the function definitions


function r = rfn(B,a,b,c,t1,t2,k,h)

        DI1 = yfnintegral(t2(1),B,a(1),b(1),c(1)) - yfnintegral(t1(1),B,a(1),b(1),c(1));
        DI2 = yfnintegral(t2(2),B,a(2),b(2),c(2) ) - yfnintegral(t1(2),B,a(2),b(2),c(2) );
        r = h*DI1/DI2 - k;
         
end

function I = yfnintegral(t,B,a,b,c)

         I1 = t.^2*(b - B*c)/2;
         I2 = -log(t)*(a*B - b*B^2/2 + c*B^3/6);
         I3 = c*t.^3/3;
         I4 = -( t*(a*B^2/2 -b*B^3/6) - a*B^3/12 )./t.^2;
         I5 = t*(a - b*B + c*B^2/2);
         I = I1 + I2 + I3 + I4 + I5;
end

[Rik]

Comment posted as a flag by Enver Can Kılıç on the previous comment:

This answer is the solution of the problem.