I'm writing a script to numerically find the integral using these three methods. Right now I can compute the right value for the integral, but I need also to print the error and the upper analytical bound of it. I have no idea how to do that in matlab. From the formulas I need to find the maximal value of f'' and f'''', but I dont know ho to do that. Can someone help me?
Thank you
Recall that the upper bound of the error for either the trapezoidal rule, Simpson's rule or midpoint rule on f(x) must be bound between the interval of x in [a,b] or else we won't be able to find this error. I have two solutions for you, depending on what has been given.
If you have a closed form solution of the integral, use the symbolic toolbox in MATLAB to first define your f(x), then use the diff command to differentiate to find f'(x). If you want the second derivative, apply another diff command to it.
Example:
syms x;
f = x^4;
df = diff(f);
d2f = diff(df);
f =
x^4
df =
4*x^3
d2f =
12*x^2
If you want to find the maximum within a certain interval, assuming a step size T between the lower bound a and b, you can try doing:
a = 1;
b = 3;
T = 0.01;
maxf = max(abs(double(subs(f,a:T:b))));
maxf =
81
What the above code does is that it finds the maximum value between [1,3] of x^4 in a step size of T = 0.01. You can replace the fourth line's first parameter f with whatever function you want. Note: The error terms require the **absolute* value, which is why the abs function is there.
If you don't have the function, but you only have a set of (x,y) pairs, you can still use the diff function on the output y values but it will implement the discrete approximation to the derivative. As such, the derivative is approximately equal to:
x_i' = x_{i+1} - x_i
This takes the (i+1) point and subtracts it with the i point. If you apply diff on an array of N points, it will return an array of N-1 points to accommodate for looking one sample ahead to find the difference. After, you can use this array and find the maximum of the derivative. If you want to find the maximum of the second derivative, simply call diff again on the first outpt of diff. Example to get you started:
x = 0:0.01:2*pi;
y = sin(x);
ymax = max(abs(diff(y)));
ymax =
0.0100
What the above does is that it finds the maximum value of the discrete approximation to the derivative to y = sin(x), as denoted by diff(y).
Hopefully one of these will work for you!
I used to teach numerical methods back when I was an instructor at my university. Take a look at these slides for a good intro to numerical integration using Simpson's Rule and the Trapezoidal Rule. There are also slides on how to code this in MATLAB (without using any built-in functions), as well as some slides on how to compute the error, and how many points you need within a fixed interval to guarantee a particular approximation error.
http://www.rnet.ryerson.ca/~rphan/mth510/f2012/Lab_10/MTH510F2012_Lab10.pdf
