1

I really need your help. I don't understand why my program works incorrectly. I need to solve this using Runge Kutta Method

y' = x*exp(pow(-x,2)) * sin(x) -2*x*y, 
y(a) = 1, a = 0,  b = 2;

Analytical answer: y(x) = exp(pow(-x,2)) * ( sin(x) - x*cos(x) + 1 )

My program:

float f(float x, float y) {
    return x*exp((-x)*(-x)) * sin(x) - 2*x*y;
}
float ya(float x) {
    return exp((-x)*(-x)) * ( sin(x) - x*cos(x) + 1 );
}

int main() {
    int n = 100;
    float a , b , h , x;
    float y[n];

    float m1, m2 , m3 , m4;
    int k;

    a = 0 , b = 2, h = (b-a) / (n-1);
    y[0] = 1;

    for (k=0; k<n-1;k++) {
        x = k * h;
        m1 = h*f(x,y[k]);
        m2 = h*f(x+h/2,y[k]+m1/2);
        m3 = h*f(x+h/2,y[k]+m2/2);
        m4 = h*f(x+h,y[k]+m3);

        y[k+1] = y[k] + (m1 + 2*m2 + 2*m3 + m4)/6;
    }

    for(k=0; k<n;k++){
        x = k * h;
        printf("x = %5.2f \t\t Numerical = %5.2f \t\t Analytical = %5.2f \n",x , y[k], ya(x));
    }
}

Result of program. Numerical and Analytical are different and I can't understand why

Evg
  • 25,259
  • 5
  • 41
  • 83
Shams
  • 11
  • 1
  • I think it is time to get out the debugger, step through all of your lines and observe the values of the variables you use. Try to spot where one variable starts to differ from your expectations. Another good approach is to put the runga-kutta (4th) order into a seperate function (e.g. like https://www.geeksforgeeks.org/runge-kutta-4th-order-method-solve-differential-equation/). – Pepijn Kramer Jun 26 '22 at 07:35
  • @PepijnKramer Please please don't refer people to that website. It's terrible. – Evg Jun 26 '22 at 08:01
  • @Evg For learning C++ you are absolutely right and the include of bits/stdc++.h and using namespace are indeed quite bad. For the runga-kutta method it was kind of fine (and in this case the code looked more like "using C++" to learn about numerical methods then about learning C++) – Pepijn Kramer Jun 26 '22 at 08:27
  • @PepijnKramer The problem is that someone with certain experience can see what's good and what's bad there, but a newbie most likely can't. If a newbie goes there to read about Runge-Kutta part, they could also easily pick up some other really bad things without realizing that. – Evg Jun 26 '22 at 08:32
  • @Evg Yes good points, and I'll take them into account – Pepijn Kramer Jun 26 '22 at 09:50

1 Answers1

2

The reason is simple: your analytical answer is wrong. Note that (-x)*(-x) is the same as x^2, not -x^2, but the derivative contains -2*x*y term. After you fix math, you'll get the expected result.

Unrelated but recommended reading: Why aren't variable-length arrays part of the C++ standard?

Evg
  • 25,259
  • 5
  • 41
  • 83