1 view (last 30 days)

Show older comments

I have a function in T which must be integrated trapezoidally. Eα is a constant. Can someone suggest how to perform such an integration in matlab?

Since this is an exponential inteagral with no unique solution, another option is to integrate by parts upto n times, ignore the nth integral and assume the terms upto n-1 as the solution. can this be done in matlab?

Thanks in advance :)

John D'Errico
on 6 Aug 2021

Edited: John D'Errico
on 6 Aug 2021

It is simply not true that there is no UNIQUE solution. There is no EXACT, analytical solution that can be written in a simple algebrac form using only addition, multiplication, etc. That is a very different thing.

It is also not true that it MUST be integrated using trapezoidal methods. For example, it is trivial to accomplish much more accurately using a tool such as integral. However, Ea MUST be given a value. Just saying that Ea is a constant is not sufficient. Any tool such as trapezoidal rule here is also meaningless if Ea is not a known constant. You cannot possibly know how to approximate the solution if you do not know the value of that constant. For example, if Ea is a relatively large number versus a very small number, then to have any degree of accuracy in the result, you will need to use very different approximations.

Doing some sort of integration by parts is silly, when you will decide to drop some terms arbitrarily. I'm sorry, but it is. I think that simply misunderstands tools like integration by parts.

But to do the computation is trivial, as long as Ea can be specified.

format long g

K = @(T,Ea) exp(Ea./(8.314*T));

intval = integral(@(T) K(T,3),270,290)

So fairly accurately, when Ea == 3, the integral is as shown. There is nothing non-unique about that solution. And there was no need to use trapezoidal rule. I could have done so. But the result would be of far lower accuracy unless we used a large number of points in the trapezoidal rule approoximation. Is the above correct? How close is it?

syms T

Ea = sym(3);

vpaintegral(exp(Ea./(sym('8.314')*T)),270,290,'reltol',1e-20)

So the result provided by integral seems accurate to as close as double precision arithmetic can report. That more than 5 significant digits is required is silly, since 8.314 is surely not known as EXACTLY that value.

What can you do if Ea is an unknown constant? The solution uses a function called ei, the exponential integral. It is a well defined function, much as is sin(x) or cos(x). Perhaps you know about sin and cos, but not about ei. That does not mean it does not exist.

syms Ea T

vpa(int(exp(Ea./(sym('8.314')*T)),270,290),5)

intEa = matlabFunction(int(exp(Ea./(sym('8.314')*T)),270,290));

And now we can evaluate that function for any value of Ea.

intEa([1 3 10 100])

Thre is no need at all to use anything fancy. Just use MATLAB.

J. Alex Lee
on 6 Aug 2021

The application still is unclear to me...

If you just want numerical integration, look at the function "integral", as implemented below.

It may be worth noting that the "exact" indefinite integral can be expressed

where Ei is the exponential integral, which seems to have a related in-built numerical algorithm in matlab: expint.

If you have crazy values of Ea, the tops and bottoms in your double sum expression may get crazy too, and then you'll have to deal with finite precision problems in the ratio. If that expression is already the result of some kind of discretization of soemthing else, may be worth going back to that and looking at the problem mathematically in more detail.

Ta = 270;

Tb = 290;

R = 8.314;

Ea = 1e+3;

num_soln = integral(@(T)intKrnl(T,R,Ea),Ta,Tb,"AbsTol",1e-16,"RelTol",1e-16);

san_soln = indefInt(Tb,R,Ea)-indefInt(Ta,R,Ea);

fplot(@(T)intKrnl(T,R,Ea),[Ta,Tb])

function out = indefInt(T,R,Ea)

Ei = -(expint(-Ea/(R*T))+1i*pi);

out = T*exp(Ea/(R*T)) - Ea/R*Ei;

end

function out = intKrnl(T,R,Ea)

out = exp(Ea./(R*T));

end

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!