Search code examples
matlabnumerical-integration

How do I make a clear table between a function and iterations, in MATLAB?

For my computing course, I am asked to determine the length of a function on the interval (0,20), by using the Trapezoidal Rule in MATLAB.

The parabolic curve is given by the following function:

y(x) = 0.1x^2 -2*x + 16

I first calculated the length through the arc length integral, by using the derivative of y: y'(x) = 0.2x - 2 I took the integral of sqrt(1+(0.2x - 2)^2). The answer should approximately be 29,5.

Now, I have the following in MATLAB:

y = @(x) sqrt(1+(0.2*x-2)^2);
x0=input('enter the value of Initial Limit x0');  *which is 0*
xn=input('enter the value of Final Limit xn'); *which is 20*
N=input('enter the total number of areas');
h=((xn-x0)/N);
area=0;
while(x0<xn)
    area=area+(h/2)*(y(x0)+y(x0+h));
    x0=x0+h;
end

fprintf('mysol for L=%.16f',area);

Now, I am asked to increase the number of used areas N in a sequence of N=200,400,800,1600... until I reach a certain accuracy.

How do I get a neat table in which I see my N vs. my length?

Solution

  • You are very close. Your code only handles the case for 1 value of N. As what Ander alluded to you above, simply wrap this code in a loop. We see that the initial value of N is 200, but you'll need to keep doubling the number until the error of your estimated area and predicted error is less than a certain amount. Be advised that 29.5 is a low precision estimate of the answer and so the value of N will be guided by this area. I would suggest increasing the precision of your value to ensure that you capture the right value of N.

    If you consult the indefinite integral on Wolfram Alpha, we get the following expression:

    enter image description here

    Therefore, we can make an anonymous function in MATLAB and find the definite integral by simply substituting 20 into this function and subtracting with 0.

    >> format long g;
>> f = @(x)((x^3 - 30*x^2 - 25*sqrt(x^2 - 20*x + 125)*asinh(2 - x/5) + 325*x - 1250) / (10*sqrt(x^2 -20*x + 125)));
>> est = f(20) - f(0)

est =

          29.5788571508919

    est contains the true value. Now we can finally modify your code. What you will have to do is when you estimate the area using the trapezoidal rule, you will need to check if the estimated value and the true value has less than some fraction of a tolerance in error. We can try something like 1e-8 where this roughly means that you will have a precision of accuracy to be about 8 digits.

    Therefore, set a variable called err to be 0.01, run this in a loop and you'll need to check if the change in error is less than this amount before you quit. You'll also want a table where we can display the number N, the estimated value, the true value and finally the error seen. We can print this off inside the loop at each iteration. Also note that you are changing the value of x0 in your estimation, so you will need to reset x0 each time in the loop:

    % Get true value of area
f = @(x)((x^3 - 30*x^2 - 25*sqrt(x^2 - 20*x + 125)*asinh(2 - x/5) + 325*x - 1250) / (10*sqrt(x^2 -20*x + 125)));
est = f(20) - f(0);

% Set the error
err = 1e-8;

% Function setup
y = @(x) sqrt(1+(0.2*x-2)^2);
x0 = 0; % Prior knowledge
xn = 20; % Prior knowledge
N = 200; % Initial value of N

% Print header of table
fprintf('%4s %11s %11s %10s\n', 'N', 'Estimated', 'True', 'Error');

while true % Start loop
    x0 = 0; % Reset each time
    h=((xn-x0)/N);
    area=0;
    while(x0<xn)
        area=area+(h/2)*(y(x0)+y(x0+h));
        x0=x0+h;
    end
    er = abs(area - est) / est; % Calculate error 

    % Print table entry       
    fprintf('%4d %.8f %.8f %.8f\n', N, area, est, er);

    % Check for convergence
    if er <= err
        break;
    end
    N = N * 2; % Double N for the next iteration
end

    When you run this, I get this table:

       N   Estimated        True      Error
 200 29.57915529 29.57885715 0.00001008
 400 29.57893169 29.57885715 0.00000252
 800 29.63483340 29.57885715 0.00189244
1600 29.60682664 29.57885715 0.00094559
3200 29.57885832 29.57885715 0.00000004
6400 29.57885744 29.57885715 0.00000001

    Therefore, this tells us that with N = 6400, this will give us the required error budget of 1e-8. This should be enough to get you started. Good luck!