Double integration with Range

Hello my friends.

I really need your help with my problem. I want to perform double integration with respect to x and y. However, x, which is the internal integration, has a range. Then, I need to plot this result.

Please note that my internal integration should be solved " numerically", I guess, as I can't find closed form solution to the internal integration, g1 in my code. Moreover, my integration are from zero to infinity.

I tried to write the code and I copied what I got so far down and I don't think this leads to the solution.

Please take a look and any suggestions will be really appreciated. Have a wonderful day!

clc
clear all;

syms a b lambda y;

x = 1:100;

    for ii = 1:length(x)

%%f1 function integration with respect to x

    f1 = @(x) lambda*exp(-a*y)*exp((b*y)/(x(ii)+1)) * exp(-lambda*x(ii));

g1(ii) = integral(f1,0,inf);

%%After finding g1 with respect to x, the desired integration is with respect to y

F(y) = (1./y).*(exp(-y) - exp(-y).*g1(ii));

lambda = sym(2);

a = sym(0.8);

b = sym(0.2);

fsub = subs(F);

Fsym(ii) = integral(fsub, y, 0, inf);

      end

      %%The plot should be with respect to x

fplot(Fsym, [0  100])

Answer 1

[Star Strider]

Try this:

syms x y
lambda = sym(2);
a = sym(0.8);
b = sym(0.2);

f1(x,y) = lambda*exp(-a*y).*exp((b*y)./(x+1)) .* exp(-lambda*x);

x = sym(1:100);

inty = double(int(f1(x,y), y, 0, Inf));

figure(1)
semilogy(x, inty)
xlabel('x')
ylabel('Integral with respect to x')
grid

[Michael Henry]

Thank you so much. Please note that I am getting the following errors:

"

Error using symengine DOUBLE cannot convert the input expression into a double array.

Error in sym/double (line 643) Xstr = mupadmex('symobj::double', S.s, 0);

Error in code (line 11) inty = double(int(f1, y, 0, Inf));

"

Can you please check it for me?

One more thing please, my request was to find the integration of function F(y) in my code above not for f1. So, can you please show me how to do that?

Thanks a lot for your time..

[Star Strider]

My pleasure.

My code works as I posted it in R2017a. I cannot tell what the problem is when you ran it.

This runs, and uses ‘F(y)’ but does not produce useful results because ‘FSym’ is uniformly infinite. I will leave that to you to sort.

The Code —

syms x y
lambda = sym(2);
a = sym(0.8);
b = sym(0.2);

f1(x,y) = lambda*exp(-a*y).*exp((b*y)./(x+1)) .* exp(-lambda*x);

x = sym(1:100);

g1 = double(int(f1(x,y), y, 0, Inf));

F(y) = (1./y).*(exp(-y) - exp(-y).*g1);
Fsym = int(F, y, 0, Inf);

figure(1)
semilogy(x, Fsym)
xlabel('x')
ylabel('Integral with respect to x')
grid

It otherwise appears to do what you want.

[Michael Henry]

I am really grateful for your kindness. It actually works now.

However, I have one thing please, the integration, g1, should be over x not y. Once we found an integration g1, we integrate over y which is Fsym in your code. When I tried to change y to x in g1 line, it gave me error that it can't integrate with respect to matrix! Can you please check it again.

[Star Strider]

My pleasure.

How does the ‘integration with respect to x’ work?

This is not clear at all from your code, or your description.

Given:

f1(x,y) = lambda*exp(-a*y).*exp((b*y)./(x+1)) .* exp(-lambda*x);

Is it:

• int(f1, x(k), Inf)
• int(f1, 0, x(k))
• Or is ‘x’ simply a parameter (constant) and not actually part of any integration except with respect to ‘y’?
• Something else?

[Michael Henry]

Thanks for your concern.

Let me explain this please: x has a range from let's say 1:100. Then, what I am trying to do is to solve the inner integral (from zero to infinity) for each point of x in that range (i.e, 1:100) and then solve the whole integral " C " which is now with respect to y. Finally, having all these solutions for every point of x, we need to plot them with respect to x itself.

I am sorry if my description still not clear. The bottom line that I want to solve double integral " C " for a range of x and then plot it with respect to this range of x.

Just a heads up, I am expecting a natural logarithmic shape for the plot.

I attached a clear image for my integration and I hope this will help more.

Please help me to solve this problem as I am really stuck for a while with it. Many thanks.

[Michael Henry]

I want to add another thing please.

If its difficult to do literally what I said in my description, can you please show me how to perform inner integration with respect to x (with y constant in it) and then perform the overall integration C with respect to y (Its basically double integral). After that, I think the plot thing will be easy as I can specify the range, like: fplot (Fsym, [0 100])? So, the range is implicitly is for x? What do you think?

Waiting your response and Thank you again.

[Star Strider]

The problem is that you want to integrate it with respect to x from 0 to +Inf. I cannot discover any way you can do that in incremental values of x, and have the result make any sense. I have been experimenting with this without success.

I am not certain that is even possible.

[Michael Henry]

What if we replace the +Inf by a very large number. Would this help?

For the results you get, Please note that this integral is not in the complete form, I just wanted to put the basic idea so it will be easier for you to deal with in your typing to MATLAB. Once you gave me the correct code, I will add the rest of the equation and I hope this will work and plot it smoothly. Thanks for your time.

[Star Strider]

I would use 1E+8 or less as the ‘very large number’.

However, I doubt that would actually work, because both your integrals are definite integrals defined as specific scalars that are the integrals from 0 to +Inf. The first remains a function of y, however the complete double-integrals are not actually functions of either x or y. The result is a scalar.

Regardless, when I evaluate them, the result is NaN.

I would welcome the insights of others here who might have ideas as to how to approach this. This is a maths problem, not the coding problem I thought it was originally.

[Michael Henry]

Can you please show me the code you have so far. I would like to read it thoroughly and see if I can do anything with it.

One more thing, is there any way to plot the integration g1 with respect to x for y, let's say, = 1?

Once again, I am really thankful for your time and effort.

[Star Strider]

My pleasure.

Not everything I wrote is relevant. My latest attempt is:

lambda = 2;
a = 0.8;
b = 0.2;
f1 = @(x,y) lambda*exp(-a*y) .* exp(b*y./(x+1)) .* exp(-lambda*x);
xv = 1:10;
for k1 = 1:length(xv)    
    f2 = @(y) integral(@(x)f1(x,y), 1E-8, xv(k1));
    If2 = f2(0);
    f3 = @(y) (exp(-y) - exp(-y) * f2(y))./y;
    If3 = f3(0);
    C = integral(f3,1E-8, 1E+4);
%     Cn(k1) = double(C);
    Cn(k1) = C;
end

I cannot figure out how to do what you requested.

[Michael Henry]

Thank you!

[Star Strider]

My pleasure!