0

Second part of Q: Then solve the integral between 0 and y of (x^2)(e^(-x^2))dx=0.1 for y using bracketing and bisection.

Here's what I have done so far:

    #include <stdio.h>
#include <math.h>
double f(double x, double y);

int main(void) {
    int i, steps;
    double a, b, y, h, m, lower, upper, x, simp, val;
    /*
     * Integrate (x^2)(e^(-x^2)) from 0 to y
     */
    steps = 20000;
    a = 0;
    b = y;
    h= (b-a)/steps;
    /*
     * now apply Simpson's rule. Note that the steps should be even.
     */
    simp = -f(a, y);
    x = a;
    for (i =0; i < steps; i += 2) {
        simp += 2.0*f(x,y)+4.0*f(x+h, y);
        x += 2*h;
    }
    simp += f(b, y);
    simp *= h/3.0;
    /*
     * print out the answer
     */
    printf("The integral from 0 to y with respect to x by Simpson's Rule is %f\n", simp);
    /*
     * Now we need to bracket and bisect to find y
     */
    lower = 0;
    /*
     * Lower bound is from turning point
     */
    upper = 100;
    /*
     *Upper value given.
     */
    while (upper - lower > 10E-10){
        m = (lower + upper)/2;
        val = f(m, y);
        if (val >=0)
            upper = m;
        if (val <=0)
            lower = m;
    }
    m = (lower + upper)/2;
    printf("The value for y is: %lf\n", m);
    return 0; 
}
double f(double x, double y) {
    return pow(x,2)*exp(pow(-x,2))-0.1;
}

Output: The integral from 0 to y with respect to x by Simpson's Rule is -0.000000
The value for y is: 0.302120

It runs but doesn't do exactly what I need it to do. I need to be able to continue working with the integral once I've used 0 and y as the limits. I can't do this. Then continue on and solve for y. It gives me a value for y but is not the same one I get if i solve using online calculators. Also, the output gave zero for the integral even when I changed the equation to be integrated to x^2. Can anyone help explain in as simple terms as possible?

Dearbhla Ryan
  • 11
  • 1
    Are you sure this is the code you use? Because it does not set `y` anywhere. – deamentiaemundi Jan 29 '17 at 00:58
  • it says:a = 0; b = y; – Dearbhla Ryan Jan 29 '17 at 01:04
  • Also there's more code if you scroll it down – Dearbhla Ryan Jan 29 '17 at 01:05
  • 2
    You cannot set `b` to `y` without setting `y` to something. – deamentiaemundi Jan 29 '17 at 01:09
  • but I don't know what y is, I need to find it, my limits are 0 and y and then I set the result equal to 0.1 to find it? Sorry if I'm not understanding. – Dearbhla Ryan Jan 29 '17 at 01:13
  • 1
    If you don't know what `y` is yet, then you simply can't do `b = y` - it doesn't make any kind of sense. What value do you expect `b` to be, after doing that? You'll need to completely rethink what you're doing, here. – Crowman Jan 29 '17 at 01:20
  • By setting `b = y;` you set `b` to some indeterminate value and the result of the whole computation is unknown. That is called "undefined behaviour". As Paul Griffiths already said: you need to attack you problem differently. – deamentiaemundi Jan 29 '17 at 01:21
  • Also, you pass `y` to your function `f()`, but that function doesn't use it, so why are you passing it? – Crowman Jan 29 '17 at 01:23
  • I see what you mean, we were given a hint that y is bigger than 0 but less than 100, but I don't know how to work with that – Dearbhla Ryan Jan 29 '17 at 01:50
  • Mmh...so you got a numerical result, say 0.422725 and are looking for the limit y (y = 2.0 for that result) numerically? – deamentiaemundi Jan 29 '17 at 01:56
  • So, I selected y=26 as beyond that point the output states that the integral is infinite and am chalking my confusion up to a misunderstanding of what is being asked of me! – Dearbhla Ryan Jan 29 '17 at 02:19
  • Online calculators also state that y≈0.745451004522108..., so there is obviously an error in the bracketing and bisection too! – Dearbhla Ryan Jan 29 '17 at 02:35
  • So your input is 0.1 and you are searching for the limit `y` (y ~ 0.74545), is that correct? Yes, you can do that with the bisection method, but your function for the bisection method is the integral itself--you compute it for several limits. Say you start at `y=26` and `y=0` you compute it for these limits, compare with your input and go on with the bisection. Might take some time (you should set a max precision to reach) but should work with this integral. If you have PARI/GP you can use `( (1-erfc(y)) * sqrt(Pi) - 2*y/exp(1)^(y^2) )/4` to check more quickly. – deamentiaemundi Jan 29 '17 at 03:18

0 Answers0