1
$\begingroup$

I'm trying to write a Mathematica function which takes an arbitrary expression as input. More concretely, I'm trying to write an integration function (utilizing the trapezoidal method), the input of which behaves like NIntegrate. For example:

fTest[x_, y_] :=  x + y;

fNew[y_] := NIntegrate[x^2*fTest[x, y], {x, 0, 10}];

Currently, I have something like:

NInt1D[f_] := 
  Module[{Num, eps, sum, a, X},
    Num = 10;
    eps = 10/Num;
    sum = 0;
    For[a = -Num, a <= +Num, a++,
      X = a*eps;
      sum += f[X]*eps;];
    Return[sum]];

Which is problematic for me, since the above NInt1d only takes functions as input, and not expressions, preventing me from defining my function in terms of the integral, as can be done with NIntegrate.

I hope this toy example makes clear my intentions. I found another post which ask this question in the context of creating a Plot function [1]. Unfortunately, the explanation was too specific, and I found it hard to understand -- I'm somewhat new to Mathematica.

[1] Passing a function as an argument of another function

$\endgroup$
  • $\begingroup$ If to fix $x$ as the integration variable, would something like this do: NInt1D[f_] := Module[{Num,eps,sum, a,g, X}, g[x_]=f; Num = 10; eps = 10/Num; sum = 0; For[a = -Num, a <= +Num, a++, X = a*eps; sum += g[X]*eps; ] Return[sum]]; ? $\endgroup$ – Andrew Dec 29 '18 at 20:27
  • $\begingroup$ You might be interested in this answer, which gives a more detailed explanation of the strategies behind implementing such a function. One important thing is that it's usually better to pass in the independent variable along with the expression, such that you're not forced to always use X for example $\endgroup$ – Lukas Lang Dec 29 '18 at 23:42
2
$\begingroup$

I studied your code for a while and it seems to me that it does not implement the trapezoid method for evaluating an integral. So even if I showed you how to make it accept both functions and expression as arguments, it would still fail to produce correct results. So, instead, I will use

f[x_] := x^2
NIntegrate[f[x], {x, -12, 12}]

1152.

for testing. Although this is simpler than the example you gave in your question, it should still serve to demonstrate a trapezoid method that meets your requirements.

I will not use For which produces ugly and inefficient code. Rather, I will use MovingAverage and Total, which provide a reasonable functoinal solution to the problem. Further, I will give two implementations, the 1st rather prolix but easy to understand; the 2nd much more concise, but it might take some study to see that it does exactly the same thing as the 1st implementation.

Method 1

Clear[nInt]
nInt[f_Symbol | f_Function] :=
   Module[{x0 = 12, nSteps, mesh, pts, eps},
     nSteps = 1000; (* should be even *)
     mesh = Subdivide[-x0, x0, nSteps];
     pts = f /@ mesh;
     eps = mesh[[nSteps/2 + 2]];
     eps Total[MovingAverage[pts, 2]] // N]
nInt[expr_, var_Symbol] :=
   Module[{x0 = 12, nSteps, mesh, pts, eps},
     nSteps = 1000; (* should be even *)
     mesh = Subdivide[-x0, x0, nSteps];
     pts = (expr /. var :> # &) /@ mesh;
     eps = mesh[[nSteps/2 + 2]];
     eps Total[MovingAverage[pts, 2]] // N]

What is ugly about the above code is that it violates the rule of "don't repeat code unnecessarily". Its main virtue is that it is easy to follow.

Method 2

Clear[nInt]
nInt[args__] :=
   Module[{x0 = 12, nSteps = 1000, argv, mesh, pts},
     argv = {args};
     mesh = Subdivide[-x0, x0, nSteps];
     pts =
       Switch[argv[[1]],
         _Symbol | _Function, argv[[1]] /@ mesh, 
         _, (argv[[1]] /. argv[[2]] :> # &) /@ mesh];
     mesh[[nSteps/2 + 2]] Total[MovingAverage[pts, 2] // N]]

Tests

nInt[f]
nInt[#^2 &] 
x =.; nInt[x^2, x]
x = 42; Block[{x}, nInt[x^2, x]]

When evaluated, all of the above expressions give 1152., the same value as Integrate gave.

$\endgroup$
  • $\begingroup$ Thank you m_goldberg. Your response is very educational -- I learned a bit more than I asked for, which I appreciate. $\endgroup$ – Paul Eugenio Jan 1 at 21:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.