0

source: https://www.ijeter.everscience.org/Manuscripts/Volume-6/Issue-4/Vol-6-issue-4-M-54.pdf

this is the code for solving the boundary value problem by the shooting method . I decided to use the formula of the secant method (or in other words, the Newton method). and decided to change the m3 search accordingly.

before: m3=m2+(((m2-m1)(b-b2))/(1.0(b2-b1)));

after: m3=m2-((m2-m1)/(b2-b1))*(b2-b);

And I get this result:

Exact value of M =-1.305508
0.500000        -1.#IND00
1.000000        -1.#IND00
1.500000        -1.#IND00
2.000000        -1.#IND00
2.500000        -1.#IND00
3.000000        -1.#IND00

that is, apparently the exact value of M is correct, but the second column does not want to count correctly. Why is this happening?I will be very grateful for your help

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
float f1(float x, float y, float z)
{
    return(z);
}
float f2(float x, float y, float z)
{
    return(x + y);
}
float shoot(float x0, float y0, float z0, float xn, float h, int p)
{
    float x, y, z, k1, k2, k3, k4, l1, l2, l3, l4, k, l, x1, y1, z1;
    x = x0;
    y = y0;
    z = z0;
    do
    {
        k1 = h * f1(x, y, z);
        l1 = h * f2(x, y, z);
        k2 = h * f1(x + h / 2.0, y + k1 / 2.0, z + l1 / 2.0);
        l2 = h * f2(x + h / 2.0, y + k1 / 2.0, z + l1 / 2.0);
        k3 = h * f1(x + h / 2.0, y + k2 / 2.0, z + l2 / 2.0);
        l3 = h * f2(x + h / 2.0, y + k2 / 2.0, z + l2 / 2.0);
        k4 = h * f1(x + h, y + k3, z + l3);
        l4 = h * f2(x + h, y + k3, z + l3);
        l = 1 / 6.0 * (l1 + 2 * l2 + 2 * l3 + l4);
        k = 1 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
        y1 = y + k;
        x1 = x + h;
        z1 = z + l;
        x = x1;
        y = y1;
        z = z1;
        if (p == 1)
        {
            printf("\n%f\t%f", x, y);
        }
    } while (x < xn);
    return(y);
}
main()
{
    float x0, y0, h, xn, yn, z0, m1, m2, m3, b, b1, b2, b3, e;
    int p = 0;
    printf("\n Enter x0,y0,xn,yn,h:");
    scanf("%f%f%f%f%f", &x0, &y0, &xn, &yn, &h);
    printf("\n Enter the trial M1:");
    scanf("%f", &m1);
    b = yn;
    z0 = m1;
    b1 = shoot(x0, y0, z0, xn, h, p = 1);
    printf("\nB1 is %f", b1);
    if (fabs(b1 - b) < 0.00005)
    {
        printf("\n The value of x and respective z are:\n");
        e = shoot(x0, y0, z0, xn, h, p = 1);
        return(0);
    }
    else
    {
        printf("\nEnter the value of M2:");
        scanf("%f", &m2);
        z0 = m2;
        b2 = shoot(x0, y0, z0, xn, h, p = 1);
        printf("\nB2 is %f", b2);
    }
    if (fabs(b2 - b) < 0.00005)
    {
        printf("\n The value of x and respective z are\n");
        e = shoot(x0, y0, z0, xn, h, p = 1);
        return(0);
    }
    else
    {
        printf("\nM2=%f\tM1=%f", m2, m1);
        m3 = m2 - ((m2 - m1) / (b2 - b1)) * (b2 - b);
        if (b1 - b2 == 0)
            exit(0);
        printf("\nExact value of M =%f", m3);
        z0 = m3;
        b3 = shoot(x0, y0, z0, xn, h, p = 0);
    }
    if (fabs(b3 - b) < 0.000005)
    {
        printf("\nThere is solution :\n");
        e = shoot(x0, y0, z0, xn, h, p = 1);
        exit(0);
    }
    do
    {
        m1 = m2;
        m2 = m3;
        b1 = b2;
        b2 = b3;
        m3 = m2 - ((m2 - m1) / (b2 - b1)) * (b2 - b);
        z0 = m3;
        b3 = shoot(x0, y0, z0, xn, h, p = 0);
    } while (fabs(b3 - b) < 0.0005);
    z0 = m3;
    e = shoot(x0, y0, z0, xn, h, p = 1);
}
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Vs_De_S
  • 155
  • 1
  • 6
  • Why did you tag this as C++? All of the code is plain C. The code could also need a couple of empty lines to divide it into paragraphs, to increase readability. And to increase understanding and maintainability of the code you should probably add some documentation (comments) as well. – Some programmer dude Dec 11 '21 at 09:15
  • Handy reading: [What do 1.#INF00, -1.#IND00 and -1.#IND mean?](https://stackoverflow.com/questions/347920/what-do-1-inf00-1-ind00-and-1-ind-mean) – user4581301 Dec 11 '21 at 09:17
  • @user4581301 yes, thanks, I have already found that this is caused by division by 0. But now I do not know how to avoid it.. – Vs_De_S Dec 11 '21 at 09:19
  • 1
    Best guess is `(k1 + 2 * k2 + 2 * k3 + k4);` summed up to 0. Stepping through with a debugger or a few diagnostic `printf`s should show you how it happened. – user4581301 Dec 11 '21 at 09:23
  • @user4581301 btw,I ran into the same problem (division by 0) with the first formula provided by the author of the code. Apparently, it is necessary to select "successful" shots here. Since in my opinion both formulas calculate the correct values. – Vs_De_S Dec 11 '21 at 09:33

1 Answers1

0

It is also wrong in the paper source, perhaps some well-meaning editor saw an apparent typo that never was.

The secant loop you want to continue while the distance to the target value is large

     do
     {
     ...
     } while (fabs(b3 - b) > 0.0005);

so that you leave this loop when the distance to the target is small enough. This then also precludes the division-by-zero.

Also, in the Runge-Kutta loop, add as first line

if(x+1.2*h>xn) h=xn-x

so that the last iteration hits the interval end exactly. This might otherwise not be the case if xn-x0 is not a multiple of h, or if floating points errors accumulate in an unfortunate pattern. This is not a concern with the test data, but could occur for h=0.1.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51