Creating the trapezoidal rule in Mathematica

[stripes]

Homework Statement

I am to design a small piece of code in Mathematica that takes in four parameters, a number z, the number of intervals n, the starting point, and the ending point. We are basically integrate a function of the form:

$f(x)=(z+1)x^{z}$

so we don't change the function at all except for the value of z.

Homework Equations

The trapezoidal rule.

The Attempt at a Solution

Here is my code so far, but my instructor is telling me that the calculated numbers/estimates are too high.

ClearAll["Global`*"]

trapRule[z_: 1, n_: 100, a_: 0, b_: 1] := 
 Module[{h, end1, end2, total},
  h = (b - a)/n;
  end1 = (z + 1)*(a^z);
  end2 = (z + 1)*(b^z);
  total = h*(end1 + end2);
  i = 1.;
  While[i <= n,
   total += 2*((z + 1)*((a + (i*h))^z));
   i++;
   ];
  (total*h)/2
  ]

I have created my own algorithms, I have looked some up online...they are all the same. I keep getting the same code above but basically it should be much more accurate for a given value of n. I can't use things like summation, I have to just use very basic commands. Again this is being done in Mathematica. Thanks in advance for any help!

 

[SteamKing]

In your handling of the endpoints 'end1' and 'end2', the value of 'total' before you enter the loop should just be the sum of 'end1' and 'end2':

total = end1 + end2

After you exit the loop, then you multiply 'total' by 'h/2'

 

[stripes]

I tried that before and handed that in, but it wasn't correct. The instructor said trapRule[4, 1000] should be the same as trapRule[4, 500, 0.5] + trapRule[4, 500, 0, 0.5], but they are not giving me the same results. He also said that I should be getting around 1.0000017 for both of those evaluations. Also, trapRule[1] should give me exactly 1 since it's just a straight line.

Perhaps I don't have the trapezoidal rule put in correctly?

 

[Bill Simpson]

Start by figuring out what step is incorrect for the simplest possible case

Trace[trapRule[0, 1, 1, 2]]

That is going to show you every step of the calculation for a single trapezoid from x==1 to x==2 of (0+1)x^0.

Then you just have to figure out what step is responsible for you getting an answer of 2 instead of 1.

Then you see if your fix also makes

Trace[trapRule[0, 2, 1, 2]]

will give you 1 instead of 1.25.

It takes a little bit of practice to understand the output of Trace[]. This
http://reference.wolfram.com/mathematica/ref/Trace.html
and
http://reference.wolfram.com/mathematica/tutorial/TracingEvaluation.html
will get you started.

 

[SteamKing]

ClearAll["Global`*"]

trapRule[z_: 1, n_: 100, a_: 0, b_: 1] := 
 Module[{h, end1, end2, total},
  h = (b - a)/n;
  end1 = (z + 1)*(a^z);
  end2 = (z + 1)*(b^z);
  total = h*(end1 + end2);   [COLOR="Red"]<-- there should be no 'h' here.[/COLOR]
  i = 1.;
  While[i <= n,
   total += 2*((z + 1)*((a + (i*h))^z));
   i++;
   ];
  (total*h)/2
  ]

What I was trying to explain in my earlier post is that you don't multiply 'total' by 'h/2' until AFTER you have accumulated all of the terms which comprise 'total'.

According to this routine, A = (h/2)*[h*(end1+end2) + all the terms between 'end1' and'end2']

 

[stripes]

SteamKing, I already removed the 'h' from that line. It's still incorrect. Trace is very confusing. I guess I'll just keep having a look at it.

 

[stripes]

I realize now that my function does the end points twice because is start at i = 1 to i = n. I have now changed it so that i = 2 and continues until i < n. But still not geting the same results for trapRule[4, 500, 0.5] + trapRule[4, 500, 0, 0.5] and trapRule[4, 1000]...

 

[stripes]

Ok, so I did the following in Python:

#!/usr/bin/env python
from __future__ import division
 
def trapRule(z, n, a, b):

    h = (b - a) / n
    s = ((z + 1)*(a**z)) + ((z + 1)*(b**z))
    for i in xrange(1, n):
        s += 2 * ((z + 1)*((a + (i*h))**z))
    return s * h / 2
 
print trapRule(4, 500, 0.5,1) + trapRule(4, 500, 0, 0.5)

and it works, so I've put the exact same code into Mathematica:

trapRule[z_: 1, n_: 100, a_: 0, b_: 1] :=
  Module[{h, total},
   
   h = (b - a)/n;
   total = ((z + 1)*(a^z)) + ((z + 1)*(b^z));
   
   For[i = 1, i <= n - 1, i++,
    total += 2*((z + 1)*((a + (i*h))^z));
    ];
   total*h/2.
   ];

trapRule[4, 500, 1/2, 1] + trapRule[4, 500, 0, 1/2]

and it appears to work, but the problem I keep running into is that Mathematica is off in its own little world again rounding the number even though i don't want it to. Python gives me:

>>> 
1.00000166667
>>>

while Mathematica gives me:

1.

When I try to use any of the following, it doesn't give me a single answer:

NumberForm[total*h/2., 10]
= 0.03125020833+0.9687514583

PaddedForm[total*h/2., 10]
= 0.03125020833+ 0.9687514583

But when I do

N[total*h/2, 10]
=1.000001667

it works very well, however, if I do:

trapRule[4, 1000]
trapRule[4, 500, (1/2)] + trapRule[4, 500, 0, (1/2)]
trapRule[4, 500, 0.5] + trapRule[4, 500, 0, 0.5]

I get:

Out[617]= 1.000001667

Out[618]= 1.000001667

Out[619]= 1.

So there doesn't seem to be any possible way to get a single answer every time. I will not get credit for the assignment if it behaves this way and this is the last assignment of my course and I cannot get credit for the course until I complete this assignment in full. It's due tonight but I can hand it in later but I use up what are called extension days (a little bit of lee way basically), and I am telling you I have been working on this same assignment for the past month to no avail. I have completed all other 18 assignments perfectly but this one, which was assigned much earlier in the semester, won't seem to cooperate. I am unsettled at the fact that Mathematica cannot evaluate such a simple and elementary algorithm.

If anyone could help me make it so that the answer is the same regardless of the format of the input parameters, I would appreciate it. Basically I want a ton of decimals every single time no exceptions, so that means if I put 1+0.5 or 1+(1/2) or 1.5 or 3/2, I want the same answer every time, no rounding at all. Thanks again.

 

[Bill Simpson]

Step 1: Get rid of this decimal point -> total*h/2.

That 2. is telling Mathematica that you are dividing by approximately 2 and you have limited precision with that 2.

Removing that gives me a result 2000003333333/2000000000000
and Mathematica by default prints floating point approximations to standard out with about 6 digits.

After you get rid of all your decimal points on input and inside your function, then if you want to see a decimal approximation of that for output then you can do
In[18]:= N[%, 40]
and get
Out[18]= 1.000001666666500000000000000000000000000
where you can see the first 7 digits are 1.000000 and Mathematica rounds that to 1.000000
and then tries to be helpful by throwing away the 000000.

You can change the PrintPrecision and a hundred other things, but for someone desperate to get their homework finished today and turned in that isn't a path to even think about starting down.

Short advice, don't use decimal points in anything, use integers and fractions and then if you really have to display a result with a decimal point use N[yourResult,nn] where nn is the number of digits you want to see.

You don't want to get mixed up using the dozen different Forms, stick with N unless you absolutely have to use something else.

Again, don't start inserting 0.5 when you can use 1/2, because 0.5 is again telling Mathematica "this is close to a half and I only have a digit or two of real accuracy in this estimate of what the number really is."

In[27]:= N[trapRule[4, 1000], 40]
N[trapRule[4, 500, (1/2)] + trapRule[4, 500, 0, (1/2)], 40]
N[trapRule[4, 500, 1/2] + trapRule[4, 500, 0, 1/2], 40]

Out[27]= 1.000001666666500000000000000000000000000
Out[28]= 1.000001666666500000000000000000000000000
Out[29]= 1.000001666666500000000000000000000000000

By the way, you can sometimes get away with putting several statements in a single cell and not separating them with semicolons, but this has been a decades long source of very hard to diagnose errors. Use Print[...];Print[...];Print[...]; or put them in separate cells. That will save you grief and confusion if you have to use this for more complex problems in the future.

If you want to put 1+0.5 and 1+1/2 and 1.5 and 3/2 to all give you the same answer every time, no rounding at all, then go get yourself another programming language. It ain't going to happen. That behavior is 25 years old and precision and accuracy and rounding and... has been a subject people have fought about for 25 years. It isn't going to change.

Ah well, you could rewrite the front-end preprocessor for Mathematica to look at all your decimal approximations and turn those into exact rationals. Perhaps surprisingly, that is possible to do. Can you learn how to do that before your 3 p.m. homework deadline? No.

 

[stripes]

Hi Bill. Thank you very much for your help. You seem very knowledgeable when it comes to Mathematica and i have learned a lot about this program in the past few months. I've used Java and Python and a few others and I understand the inherent problems that many beginners like myself will run into like precision and what not, but it's not something I understand completely. Dividing by 2. should be no different than dividing by 2. If I want it to spit out this many decimals it should always do that no matter what. The programmer should always have complete control over the code they are using. But alas, that is not the case here. And that is not the case with most languages.

I'm also frustrated that there is basically no undo button in Mathematica. That and there are an insane and ridiculous amount of brackets that can be confusing. Put these two things together and it's a disaster.

Anyways, so you're suggesting I leave it as

trapRule[z_: 1, n_: 100, a_: 0, b_: 1] :=
  Module[{h, total},
   
   h = (b - a)/n;
   total = ((z + 1)*(a^z)) + ((z + 1)*(b^z));
   
   For[i = 1, i <= n - 1, i++,
    total += 2*((z + 1)*((a + (i*h))^z));
    ];
   total*h/2
   ];

so it returns rational numbers for the examples I talked about? And then just use N[expr,m] for displaying decimals?

 

[stripes]

When I first had the assignment marked the instructor said ""trapRule[4, 1000]" is not returning the same value as "trapRule[4, 500, 0.5] + trapRule[4, 500, 0, 0.5]". In fact you should be getting something close to 1.0000017 for these."

What is annoying is that I have not yet been able to give the same result for these two cases no matter what I try, and i will surely not get credit for this assignment if I can't do this...

Edit: hopefully I can suggest to the instructor that i use 1/2 instead of 0.5

 

[Bill Simpson]

Dividing by 2. should be no different than dividing by 2. If I want it to spit out this many decimals it should always do that no matter what.

That is your opinion about how a programming language should behave. There are many different programming languages with many different models behind them. Really learning half a dozen very very different languages with completely different underlying mental and calculation models can give you some perspective, but you don't have time for that now.

The programmer should always have complete control over the code they are using. But alas, that is not the case here. And that is not the case with most languages.

That has not been the case for decades. Take the gigabytes of code, good and terrible, making up your operating system, decompile that, reverse engineer that, and get t a point where you have complete control. That isn't ever going to happen.

I'm also frustrated that there is basically no undo button in Mathematica. That and there are an insane and ridiculous amount of brackets that can be confusing. Put these two things together and it's a disaster.

What I think you really really want is the "do what I mean" button that figures out what you are thinking and makes it do things the way you want. Write your own language.

This is probably a very bad thing to hand you and is opening up the door for even more problems that you won't expect.

In[1]:= trapRule[z_: 1, n_: 100, a_: 0, b_: 1] :=
Module[{rationala, rationalb, h, total},
rationala = Rationalize[a];
rationalb = Rationalize;
h = (rationalb - rationala)/n;
total = ((z + 1)*(rationala^z)) + ((z + 1)*(rationalb^z));
For[i = 1, i <= n - 1, i++,
total += 2*((z + 1)*((rationala + (i*h))^z));];
N[total*h/2, 40]];
Print[trapRule[4, 1000]];
Print[trapRule[4, 500, 0.5] + trapRule[4, 500, 0, 0.5]];

From In[2]= 1.000001666666500000000000000000000000000

From In[3]= 1.000001666666500000000000000000000000000

Then you have to decide what happens when Rationalize doesn't give you the nice exact value you think it should. Or when z isn't an integer or rational. Or when n isn't an integer... What happens when you have a number like 1/3 that goes on forever in decimal or a number that has different problems in base 10 than it has in base 2 (that is used for floating point calculations). The list goes on and on. What happens if the result needs more than 40 digits? What happens if he doesn't like 40 digits and thinks there should be less?

What should really be done with your code, I can't decide what path to take.

 

[stripes]

Of course it is my opinion and this is why I will never go into programming. I would like to have it as a skill but math is my calling (I think). Of course math is the basis for computer science..but anyways.

I don't think there's a soul on this planet who would say no to a do-exactly-what-i-am-thinking-button, but when I say I wish there was an undo button, I mean exactly that. Often times I will cut the wrong part out because the mouse pad on my laptop is jittery and then I'll do one more action and then I cannot undo past that one action. Meaning I can't get back anything I had before the last action. Sometimes I can't even undo the most recent action. This means I have to close and start over from where I saved last.

That has easily been the most frustrating thing about Mathematica. Fortunately this copy was free, but had I paid for it, I would be incredibly mad, just because of the undo button alone. It has been a powerful tool but it can act very strange sometimes.

I appreciate your help but your code won't be necessary. I handed in what I had because he said my code had looked alright so far. It does what it needs to do and thanks for your help.

I'm sure I'll be back if there are more problems but not to copy and paste your code, but to figure out what is wrong. Thanks again.

 

[stripes]

So I just got the mark back and I got credit for the assignment. I just ended up doing the calculations without any decimals or N[] or anything and then did N[%, 10] at the end. Thanks again Bill.

 

[Bill Simpson]

In mathematics, as opposed to Mathematica, there is an agreement that reader and writer both have "mathematical maturity" which roughly means "I'll give you what I think are enough hints and you figure it all out from the beginning for yourself and overlook any small errors or inconsistencies and think that those are not really important and just an annoyance."

Mathematica has no "mathematical maturity", in fact I don't believe anyone has ever written a piece of software which has that. I remember reading a few years ago about the various formal proof verification software tools. Those will take a sketch of a proof and declare they have been able to verify it or not. The claim was that after becoming highly skilled at using one of these a person would be able to verify a single page of an undergraduate math text in about a day of intense effort. It takes that much time and effort to provide the level of detail and information and correctness needed for the tool to be able to verify the proof. And those are probably the best tools that anyone knows how to construct today.