Non-integrable multiple integrals for Mathematica

[Anna Kaladze]

“Non-integrable” multiple integrals for Mathematica

Dear all,

I have been trying to crack one problem in Mathematica, but I keep getting a wrong answer probably because I have something either fundamentally wrong analytically or code wise. OK, here is the problem.

Suppose we have to evaluate the following:

Integral from 0 to t, dm [of Integral from m to t, dd [of Integral from d to b of (y-d), dy],

and suppose we know that b = 3. Clearly, it is straightforward to accomplish this task (please see the first 3 steps in the code attached (Example2.nb files) and the “final” answer). (I deliberately proceeded in simple steps to keep a better track of the individual answers).

Suppose also I know that t is a time variable which varies from, say, 0 to 2. Then I can easily take various values of t from that interval, plug them into the “final” answer and plot if I need to. (Step size 0.01 would be precise enough for my purpose).

However, suppose we want to plot the “final” for "all" possible values of t, without being able to analytically integrate the first inner integral, i.e. “Integral from d to b of (y-d), dy”. Let us just pretend that the reason is that this thing is “non-integrable”.

I have performed a sequence of simple steps to get to the “finalnew”, but you can see that this is a wrong array, because if I set t=0.01, for example, in “final” I would get 0.000224001, which is not what the 2nd element in “finalnew” array is.

I would really appreciate if you could help. This is for my research, where the primary problem is that in place of (y-d) I have something pretty-complicated, which is “non-integrable”. I can numerically integrate it given that d was a dummy for t which goes from 0 to 2 as you recall. But then the question is, given 2 remaining integrals, how to build a final array of numbers to plot. All I need is the final array of right numbers. I am sure there is a correct and perhaps an easier way to accomplish this task, but the question is what it is.

Thanks a lot for your time.

Anna.
-----------------------------------------------

 

[Bill Simpson]

Would it be possible for you to make the second half of your notebook significantly less complicated?

I have made a serious effort to go through your notebook three times, each time starting from scratch trying to decompile it back into something much simpler that I can understand and thus figure out where the error is. And I have failed.

You do not have any outright syntax errors, which is better than many questions posed here, but you have integrals of lists of constants that you transpose into matricies only then to turn back into vectors by applying parts of the row to the rest of the row.

Your example in the first half, where you can use Integrate, simplifies down to Integrate[y - d, {m, 0, t}, {d, m, t}, {y, d, b}]

For the second half it is clear you have spent a lot of work trying to find a way to force the available Mathematica tools into a form that will try to calculate what you want. But is there any more conventional Mathematica construct you could use instead of repeatedly using Apply on the Transpose on the Integrate. I'm thinking this must represent something simpler, perhaps possible to create with something as simple as a Table[NIntegrate[something],{range},{range}].

If not then perhaps a simple English description of what you are trying to accomplish with each step might make it feasible to find a simpler and more easily understood way of doing this that would have as a side effect also being correct.

Or as a last resort, could you just manually calculate and show exactly what the values in the first half dozen or dozen entries of finalnew should be and precisely how you calculate the corresponding different values in final?

At one time I had almost reached the conclusion that part of the problem might be that your first NIntegrate always uses b as the upper limit of integration, rather than integrating only over an interval of width step, but after another attempt from scratch to decompile what you are really trying to accomplish I am not even convinced that is correct.

 

[Anna Kaladze]

Hi Bill,

Thanks for your reply.

I have, BTW, just figured out what my error was. I failed to take into the account the fact that the first step should represent not just a single number, but the multiple areas under the curve. My approach was mathematically wrong.
Unfortunately, that confused me even more, frankly, i.e., now I am unsure how to crack the problem.
I have attached a pdf file, which explains the problem, and an nb file, where I have just assumed that F(y,d) is y-d for the sake of illustration. All I need is get to the answer in the form of a “finalarray” as shown, but without using Integrate (assuming that y-d cannot be Integrated directly).
I have tried using interpolating function, etc. but my poor programming skills result in being stuck.
I appreciate your help.

Anna.
P.S. This is NOT a homework assignment, but this just would greatly simplify my research as I deal with this type of problem a lot.
P.P.S. I have deliberately didn't insert into an nb file an original function F(y,d) since I know for sure it is "non-integrable", and it would occupy like half a page. If I know the procedure for a simple y-d, then I can apply it for other functions.

 

[Bill Simpson]

You ask for the final answer in the form of an array. Is that array containing the integral for various values of the unknown t? The unknown m? (I'm guessing not) I'm trying to clearly identify exactly what the index of that array represents, understand exactly the form of the "prism" that you are integrating over and clearly see how your function F relates to the individual dVolume bits making up the prism.

Not showing pages of meaningless detail for your F(y,d) is perfectly fine, better than confusing the issues by showing that. If, for the moment, you just do not define f[y,t] you could insert that into your Mathematica code and not try to actually evaluate the result, but just show what you are hoping to calculate.

NIntegrate, like the other Nfunctions... NSolve and NDSolve, is strict in its' demands that everything evaluate to numeric constants. You cannot expect NIntegrate to evaluate any expression that contains free variables.

I assume, and please confirm this, that f[y,d] is not independent of y and d, otherwise why would you have parameters y and d, but is independent of m and t. Thus, in light of the previous paragraph, I believe you cannot use NIntegrate on f[y,d] for a single value of y and d and then use that result across all values of y and d. But if my assumption is correct then you can NIntegrate your innermost f[y,d], that would result in a scalar constant and the two outer integrals appear to be over a triangular area, so you multiply that scalar times the area of the triangle and you would be done.

But your second level integration involves dd and I can't quite see yet how that is going to interact with your innermost integration.

Can you explain some of the things I still do not understand?

 

[Anna Kaladze]

Hi Bill,
Thanks for your response.

1. Yes, you are right, I need a final array of numbers for various values t. From the pdf file if we assume that F(y,d) is y-d, then (1) results in 9t^2/4 –t^3 + t^4/8. I do not need to “see” what the final explicit answer is in terms of a mathematical function of t, but I do need to know the “approximate” values of the function for the values of t from 0 to 2 in step size 0.01.
2. Yes, f(y,d) depends only on y and d, and not on m, t or anyhting else.
3. By assumption, b is given and d takes trhe values from 0 to t(=2) and F(y,d) is such that it is possible to aplly NIntegrate.

I think my problem is the coding issue. That is, given the above assumptions, how to code the problem which would construct that final array of 201 elements, without applying Integrate command.
Thanks.
Anna.

 

[Bill Simpson]

Just out of curiosity, exactly what does the following give you when you evaluate it?

f[y_, d_] := yourlargefunctiondefinition;
Table[NIntegrate[f[y, d], {m, 0., t}, {d, m, t}, {y, d, 3.}], {t, 0., 2., .25}]

 

[Anna Kaladze]

Dear Bill,
Thank you very much for your reply. Your method is indeed very efficient and clever, and I have got the results, all being real numbers, with no hassle in the computing NIntegrate process...I appreciate the time you have invested into this question.
I have tried to extend your method to a “slightly” different problem, which I also often face, but perhaps because I making a logical/procedural error, I get stuck.
Can you please see the attached pdf file? Essentially, the problem “looks” similar to the previous one, but this time the inner integral is an exponent of a double integral. All other assumptions remain the same. I am attaching the nb file too.
I tried various combinations, and even Log substitution for an inner integral, but nothing helps...Perhaps I need to develop a better understading of how Mathematica handles multiple integrals, but so far I didn't come across to a problem very similar to this one ...
Thanks a lot,
Anna.

 

[Bill Simpson]

I'm glad we were able to find such a simple solution to your initial problem. Just to be cautious, you should carefully check the numerical results of my NIntegrate expression to make sure that it is giving you right answers.

 

[Dale]

NIntegrate is definitely the way to go, and it is much better to do it as a single command than to split it up into three integrals. The numerical routines can control the overall accuracy and precision if you evaluate it as a single NIntegrate. Be sure to read the documentation on NIntegrate so that you understand what it is doing.

 

[Anna Kaladze]

Hi Bill,
Thanks a lot for your reply.
Just to make sure that I understood you correctly, you are asking me to verify if (1) (please see the attachment) is true or not? I will be working on it to see if I can verify it.
Please let me know.
Thanks a lot.
Anna.

 

[Bill Simpson]

I have made an error in what I was trying to accomplish.

 

[Anna Kaladze]

Now as an alternate formulation multiply each dy*dd*dm by E^m, integrate over dy and dd to get E^m times the area and finally integrate over dm to get your weighted volume.

Hi Bill,
Sorry for a silly question, - I have been trying various techniques to get this problem into a simpler looking triple integral according to your suggestions, - so I might be a bit slow ... When you say "multiply each dy*dd*dm by E^m", do you mean that you "insert" E^m into each of the 3 integrands? And then you use the integration by parts for every innermost integral?
Can you please clarify?
Thanks a lot.
Anna.

 

[Anna Kaladze]

Bill,
Correct me if I am wrong, but it appears to me that when we try to make E(of double integral) become of the form of a double integral (of some function), then essentially weare solving an integral equation. If so, I fear that this might not be an easy task...But I might be wrong ...
Thanks a lot.
Anna.

 

[Bill Simpson]

Ah, no. I am sorry but I believe I have led you down a path that does not work. I think I see now that the idea is incorrect and cannot be repaired.

I am very sorry for my error and taking up your time with that.

So I think you are back to your Integrate[E^Integrate[]] problem. If you can try to show exactly how you are having a difficulty evaluating that then perhaps we can find a solution for you.

 

[Anna Kaladze]

Ah, no. I am sorry but I believe I have led you down a path that does not work. I think I see now that the idea is incorrect and cannot be repaired.

I am very sorry for my error and taking up your time with that.

So I think you are back to your Integrate[E^Integrate[]] problem. If you can try to show exactly how you are having a difficulty evaluating that then perhaps we can find a solution for you.

Dear Bill,
Please do not apologize as I really appreciate your time and very useful hints.
I will try to formulate a solution, and see if it can get us somewhere...I will be in touch soon. My only concern now is as follows: since the problem is back to the original, i.e. cannot be represented as a triple integral, do you think I must use some kind of nested looping to get to the answer, or you believe a proper representation of NIntegrate[etc.etc...] can still work?
Many thanks.
Anna.

 

[Dale]

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.

 

[Anna Kaladze]

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.

Dear DaleSpam,
Thanks you very much for your response.
In my problem (please see the attachment) t is time variable, which is assumed to go from 0 to 2 (step size in 0.01) is sufficient for the accuracy.
That is the final answer is an array of 201 elements. I wonder if your method can be used where t goes from 0 to 2 in step 0.01. Maybe combining it with a Table command?
I am playing with it now to see how far I can get...
Thanks a lot.
Anna.

 

[Anna Kaladze]

P.S. Hi DaleSpam,
I have tried this Table[NIntegrate[
Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f ->
Function[{y, d}, y], b -> 3}, {t, 0, 2, 1/4}], but I fear I might be wrong...

 

[Dale]

That ran on my computer.

One other thing to think about. If your function is fairly smooth then you may have a computational advantage by making an interpolating function for the argument to the Exp function. However, I am not certain of that.

 

[Anna Kaladze]

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.

Hi DaleSpam,
Sorry for a silly question, but what does {f -> Function[{y, d}, y] exactly accomplish in the ab ove code? It appears that without the task can be accomplished too, but maybe I misunderstand something...Can you please also advice where I can read more about that syntax?
Many thanks.
Anna.

 

[Anna Kaladze]

Dear DaleSpam and Bill,
Thanks a lot for all your responses as they have been extremely helpful and made me learn a lot. I appreciate both of yours time.

To make sure that I have a proper understanding, I started twisting the problem a little bit (and also because exp problem is triky to symbolically integrate), and it seems to me that the method you have offered works.

If you see the pdf file and nb file, and it looks like it works fine and it given me an error (please see “diff” value less than 1/1000000).

One thing I am worried about is the property of F being well-behaved, as DaleSpam has indicated that the method would work in such a case. I wonder if the method would not for more complex functions...
Also, can you please help me to understand

1) why so I get those “\[NoBreak]y\[NoBreak] = \[NoBreak]d\[NoBreak] is not a valid limit of integration. ” errors and if they pose any “danger” really? I also wonder why despite those errors the right answer is reached anyway.

2) Further to my very last post, as I have deleted f-> function ... syntax, I wonder if this is OK to do. Maybe yes for this particular one, but not the one where had this exponent of double integral?

Thanks a lot.
Anna.

 

[Dale]

Sorry for a silly question, but what does {f -> Function[{y, d}, y] exactly accomplish in the ab ove code?

That is just where the function f(y,d) is defined. For testing purposes I just had f(y,d)=y.

Can you please also advice where I can read more about that syntax?

F1, the online help.

 

[Dale]

1) why so I get those “\[NoBreak]y\[NoBreak] = \[NoBreak]d\[NoBreak] is not a valid limit of integration. ” errors and if they pose any “danger” really? I also wonder why despite those errors the right answer is reached anyway.

Those are due to the order of evaluation. First it tries to evaluate the NIntegrate, it finds that there are still some symbolic limits, and throws the error. Then it makes the substitution, re-evaluates NIntegrate, finds that there are no symbolic limits and succeeds.

 

[Anna Kaladze]

Hi DaleSpam,

It seems the method you have suggested is really flexible and efficient and can be applied to a number of variations succesfully!
Unfortunately, as those variations get messy enough, it becomes impossible for me to re-check whether the answer is exactly the same as “Integrate” commands would imply (especially with Exp).
I have just wrote a small code (please see the nb attachment) for the problem in the pdf file. However, I fear that maybe I have made a mistake due to a potential misunderstanding of your idea. Can you please have a quick look to see if what I wrote is indeed a correct representation of the problems solution, and does not go against your idea? The code obviously runs well, but I fear that I have made an error.
I appreciate your time.
Anna.

 

[Dale]

It looks pretty good. Here are a couple of small suggestions.

Because you are repeating that inner-most integral so many times it is actually quite a bit more efficient to evaluate it once, construct an interpolating function and then just evaluate that interpolating function in the other integrals. It calculates about 26 times faster with little loss in accuracy.

 

[Anna Kaladze]

It looks pretty good. Here are a couple of small suggestions.

Because you are repeating that inner-most integral so many times it is actually quite a bit more efficient to evaluate it once, construct an interpolating function and then just evaluate that interpolating function in the other integrals. It calculates about 26 times faster with little loss in accuracy.

Hi DaleSpam,
Thanks a lot for your reply and suggestions.
Can I ask for some clarifications?
When you say the innermost integral, do you refer to the one where the limits are d and b?
And then you suggest evaluating it once and then interpolating. Can you please explain how would I evaluate it? Because the limits of integrations, d, m and t are all variable, and I am not sure how to do it. If I know, I can as you suggest construct an interpolating function, but I am not sure. Can you please help?
Thanks a lot.
Anna.

 

[Anna Kaladze]

Hi DaleSpam,
I might have figured out what you have suggested. Can you please have a look at the attached nb file and let me know if I understood you correctly. Or if there is even more efficient way to run this procedure for my problem? I wonder if what I wrote makes things faster actually...
Thanks.
Anna.

 

[Dale]

For function interpolation you will be better off using the example I showed in post 25 than the one you used in post 27. The built-in function will look at the function it is interpolating and make sure to adapt the sampling density if the function is changing rapidly.

 

[Anna Kaladze]

For function interpolation you will be better off using the example I showed in post 25 than the one you used in post 27. The built-in function will look at the function it is interpolating and make sure to adapt the sampling density if the function is changing rapidly.

Thanks a lot for your help, I have not noticed your file in post 25.

I have used now this approach, and indeed it is better.

For a different function and a slightly modified ranges of the integration limits, I after I interpolatefunction, the following message:

FunctionInterpolation::ncvb: FunctionInterpolation failed to meet the prescribed accuracy and precision goals after 6 recursive bisections near j={0.}. Continuing to refine elsewhere.

I wonder if I should be ignoring this message or if it can cause any harm.

Thanks a lot.

Anna.

 

[Bill Simpson]

I tend to be cautious and I think it would be very nice if someone would contribute the perhaps small changes needed to get rid of all the error and warning messages. Sometimes those are there for good reasons and I often learn something when I find out why.

 

[Dale]

You probably have a division by 0 somewhere. Look at your variable j and see what happens near 0. Possible solutions are to interpolate the inverse of the original equation, or to simply avoid 0.