Fourth order Runge–Kutta in C# - two differential equations

[Itosu]

Hello!

I am trying to solve two differential, non-linear equations in C# using fourth-order Runge-Kutha method.


I implemented them in C#:

 private const int ALFA = 7600;
        private const int MU = 7200;
        private int BETA = 1430;//int.Parse(HttpContext.Current.Session["BETA"].ToString());
        private const int LAMBDA = 2470;
        private int V = 139000;//int.Parse(HttpContext.Current.Session["V"].ToString());
        private const double FI = 0.51;
        private const double TETA = 2.5;
        private const int QL = 8400;
        private const double H = 0.01;

        public double x = 0.81;
        public double y = 0.055;
        // ;
        public ArrayList yList;
        public ArrayList xList;
    

        public void Start()
        {
            xList = new ArrayList();
            yList = new ArrayList();
            double t = 0;
            
            //   xList.
            // y = Y0;
            for (double krok = 0; krok <= 500; krok++)
            {
                t = krok / 100; 
                xList.Add(x);
                yList.Add(y);
                x = rungeDX(x, t);
                y = rungeDY(y, t);
                //x = X0 + (i * H);
                //y = runge(x, y);

            }

        }

        public double doplyw(double t)
        {
            if (t >= 0.5 && t < 0.75)
                return 100000;
            else return 0;
        }


        private double dx(double x, double y, double t)
        {
            double result;
            if (x > TETA)
                result = (doplyw(t) + QL - LAMBDA * x - V * x * y - MU * (x - TETA)) / 15000;
            else
                result = (doplyw(t) + QL - LAMBDA * x - V * x * y) / 15000;
            return result;
        }

        private double dy(double x, double y)
        {
            double result;
            if (x > FI)
                result = (-ALFA * y + BETA * (x - FI)) / 15000;
            else
                result = (-ALFA * y) / 15000;
            return result;
        }

        double rungeDX(double x, double t)
        {
            double K1 = dx(this.x, y, t);
            double K2 = dx((this.x + 0.5 * H), (y + 0.5 * K1), t);
            double K3 = dx((this.x + 0.5 * H), (y + 0.5 * K2), t);
            double K4 = dx((this.x + H), (y + K3), t);
            double rest = (K1 + 2 * K2 + 2 * K3 + K4) * H;
            double runge = this.x + rest;
            
           
            return runge;
        }

        double rungeDY(double x, double t)
        {
            double K1 = dy(x, y);
            double K2 = dy((x + 1 / 2 * H), (y + 1 / 2 * K1));
            double K3 = dy((x + 1 / 2 * H), (y + 1 / 2 * K2));
            double K4 = dy((x + H), (y + K3));
            double rest = (K1 + 2 * K2 + 2 * K3 + K4) * H;
            double runge = y + rest;
            return runge;
        }

    }

However the results are not that that should be (in comparison to Simulink model).
Do you have any ideas?

Thank you in advance!

 

[I like Serena]

Welcome to PF, Itosu!

You have mixed up to the function of x and t in your formulas.

Btw, it's a bit weird to have an argument x to a function if you're not going to use it, but use this.x instead.

Anyway, your code

        double rungeDX(double x, double t)
        {
            double K1 = dx(this.x, y, t);
            double K2 = dx((this.x + 0.5 * H), (y + 0.5 * K1), t);

should be something like:

        double rungeDX(double x, double t)
        {
            double K1x = dx(this.x, y, t);
            double K1y = dy(this.x, y, t);
            double K2x = dx((this.x  + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);

 

[Itosu]

Newer version below:

 

[Itosu]

The newest version is like this:

        public void rungeDX(double x, double y, double t)
        {
            double K1x = dx(x, y, t);
            double K1y = dy(x, y);

            double K2x = dx((x + 0.5 * H), (y + 0.5 * K1x * H), t + 0.5 * H);
            double K2y = dy((x + 0.5 * H), (y + 0.5 * K1y*H));

            double K3x = dx((x + 0.5 * H), (y + 0.5 * K2x * H), t + 0.5 * H);
            double K3y = dy((x + 0.5 * H), (y + 0.5 * K2y * H));

            double K4x = dx((x + H), (y + K3x * H), t + 0.5 * H);
            double K4y = dy((x + H), (y + K3y * H));            double restX = (K1x + 2 * K2x + 2 * K3x + K4x) * H / 6;
            double restY = (K1y + 2 * K2y + 2 * K3y + K4y) * H / 6;
            this.x = x + restX;
            this.y = y + restY;        }

and looks like the one in Simulink :)
any comments? what about "t"?

 

[I like Serena]

The newest version is like this:

and looks like the one in Simulink :)
any comments? what about "t"?

Yes, I have some comments. :)
It should be like this:

        public void rungeDX(double x, double y, double t)
        {
            double K1x = dx(x, y, t);
            double K1y = dy(x, y);

            double K2x = dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
            double K2y = dy((x + 0.5 * K1x), (y + 0.5 * K1y));

            double K3x = dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
            double K3y = dy((x + 0.5 * K2x), (y + 0.5 * K2y));

            double K4x = dx((x + K3x), (y + K3y), t + H);
            double K4y = dy((x + K3x), (y + K3y));


            double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
            double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;



            this.x = x + restX;
            this.y = y + restY;


        }

Cheers!

 

[Itosu]

Thank you for your response, but unfortunately I do not understand why :)
My solution gives me exactly the same graphs that a simulation in Simulink, and your's does not.
However, you are the master, not me, so I am afraid that you are correct, but please explain :)

 

[I like Serena]

Thank you for your response, but unfortunately I do not understand why :)
My solution gives me exactly the same graphs that a simulation in Simulink, and your's does not.
However, you are the master, not me, so I am afraid that you are correct, but please explain :)

I checked the http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods" and I see we need another modification.
It should be:

        public void rungeDX(double x, double y, double t)
        {
            double K1x = H * dx(x, y, t);
            double K1y = H * dy(x, y);

            double K2x = H * dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
            double K2y = H * dy((x + 0.5 * K1x), (y + 0.5 * K1y));

            double K3x = H * dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
            double K3y = H * dy((x + 0.5 * K2x), (y + 0.5 * K2y));

            double K4x = H * dx((x + K3x), (y + K3y), t + H);
            double K4y = H * dy((x + K3x), (y + K3y));


            double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
            double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;



            this.x = x + restX;
            this.y = y + restY;
        }

You can see the regular form on wikipedia, where they define the problem as:
y' = f(t, y)
and they use a time step h.

What you need to be aware of, is that this y is actually the vector (x,y).
So everywhere where you read "y" in the wikipedia article, you need to replace it by "x, y".
The intermediary constants k1, k2, k3, k4 are also vectors, meaning they have an x and an y component.
And instead of a function f(t, y) you have 2 functions dx(x, y, t) and dy(x, y), which together form a vector again.

Furthermore in your program you have reversed the order of the position vector and the time.

 

[Itosu]

You are totally correct. Thank you!
Do you know of any possibilities to calculate an error of such solution?

 

[I like Serena]

You're welcome.

The typical method to estimate the error is the following.

1. Calculate the step with stepsize H yielding a new point (x1, y1)
2. Recalculate the step with stepsize H/2, and then calculate another step with H/2, yielding a point (x2, y2).
3. The difference between these 2 points is the estimation of the error.

The resulting error is accurate up to O(h⁵).

This method is typically used for the adaptive step size Runge Kutta method.
In this method the step size is dynamically adjusted (halving it or doubling it) to keep the error within a preset bound.

 

[Itosu]

Hi :)
Thank you again!

And one more question...
I would like to calculate Area Under the Curve that I printed after having resolved the equations. I know I can do this with rectangle method, adding the areas of each rectangle, however I have a premonition that it has huge error then. What do you think?
Is there another possibility to calculate AUC?

and.. thank you in advance!

 

[I like Serena]

Rectangle method would work if you take your step size small enough.
To increase the accuracy you need to weigh your function evaluations differently.

On wikipedia you can find a list of methods:
http://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas"

You would want to use the composite rule, meaning you subdivide your interval in many small intervals, and combine that with one of the schemes listed.If f₀, f₁, f₂, ..., f_n are your function values with a step size H between them.
You have for instance the options:

Trapezoid rule (error is O(H³) which is less accurate than Runge-Kutta).
n * H * (f₀/2 + f₁ + f₂ + f₃ + ... + f_(n-1) + f_n/2)

Simpson's rule with n even (error is O(H⁵) which is more accurate than Runge-Kutta):
(n * H / 6) * (f₀ + 4 f₁ + 2 f₂ + 4 f₃ + 2 f₄ + ... + 2 f_(n-1) + f_n)Of course you can make it as complex as you'd like, but I guess first you would need to find out what problems you run into exactly if any.
If you have discontinuities or wild fluctuations in your derivatives you may need to take special care.
As it is I, believe you have a discontinuity, so you may consider doing the integration in such a way that you first integrate up to the discontinuity if that's feasible, and then add a second integration starting from the discontinuity.

If necessary, there are adaptive integration methods that may help, but if you're programming this yourself, I would only go there if you have need.
As an alternative, you may consider using math packages.
There are bound to be a few available for C#.NET.

 

[AlephZero]

Hi :)
Thank you again!

And one more question...
I would like to calculate Area Under the Curve that I printed after having resolved the equations. I know I can do this with rectangle method, adding the areas of each rectangle, however I have a premonition that it has huge error then. What do you think?
Is there another possibility to calculate AUC?

Erm... "calculating Area Under the Curve" is the same as "integrating".

You just got your code for integration working, so why do you want to use a different method for the next part of the problem?

 

[Itosu]

What do you mean? Can you explain?

 

[I like Serena]

You did not really explain what you need the area under "the curve" for.
Depending on what you want exactly, it's possible to let Runge-Kutta do the work for you.

What you have is functions dx and dy that define a vector field.
You're using RK to find a curve with a gradient that matches these functions.
Effectively you are integrating the functions dx and dy.

If I understand you correctly, you want to know the area under the curve you just found.
This means that you are integrating the function again.
You could set this up as a new second order differential equation, using RK to find the result, instead of the first order differential equation that you are solving now.

It depends on what your results are for, to decide what's best and easiest in your case.
If you really only want the curve, and only want to display the area under the curve as an informative extra, I would integrate as I described in my previous post.
If you actually want the area under the curve, you can set up a different calculation and let RK find that for you.

 

[Itosu]

My results are for calculations - I need to know the exact value od AUC (not for printing, just for forward calculations).
Do you mean that if I use the same (4th order) RK solver I will get AUC? Could you provide more details?

Thank you very much!

 

[I like Serena]

Okay, I've just worked out the math, so here goes.If I've understood everything correctly, what you would need to do is introduce two extra variables in your class:
Y (to keep track of the surface) and
t (to keep track of time).
Furthermore you would have to rewrite the rungeDx step function (shall we call it the rungeStep function?), and adjust the Start function.
And that's all!

You may not understand everything I write below, but could you at least answer the questions set in bold?

@AlephZero: Feel free to comment.
First let's rename your functions dx and dy to fdx and fdy to avoid ambiguity.
So you would have functions fdx(x, y, t) and fdy(x, y) that calculate the time derivatives of respectively x and y at some point (x, y) at time t.

Your current problem is then defined by de differential equation (DE):

dx/dt = fdx(x, y, t)
dy/dt = fdy(x, y)​

Btw, the only dependency on t, is apparently an impulse function for 0.5 <= t < 0.75 that is added to fdx(x, y, 0).

Until now you used RK to find a solution (x(t), y(t)) given a starting point (x₀, y₀, t=0).

If I understand correctly, what you're actually looking for is the function Y(x) which I define to be the area under the curve y(x) between x=x₀ and some value of x.
Can you confirm this?
To find Y(x) with RK we need to set up a new DE which is:

dY/dx = gdY(x, Y)​

So the problem becomes finding this function gdY, and also take care of the impulse dependency on t.

We can do this by hauling y back in again, which is the derivative of Y(x), and also getting t back:

dY/dx = y
dy/dx = gdy(x, y, Y, t)
dt/dx = gdt(x, y, Y, t)​

Differentiation of y(x) gives:
dy/dx = (dy/dt) / (dx/dt) = fdy(x, y) / fdx(x, y, t)

Your problem would be solved by feeding the following system to RK:

dY/dx = y
dy/dx = fdy(x, y) / fdx(x, y, t)
dt/dx = 1 / fdx(x, y, t)​

There are a few boundary conditions that have to be satisfied though.
1. t must be a well defined function of x.
2. y must be a well defined function of x
3. dx/dt must never be zero.

So here I'm getting back to you.
Can I assume that x strictly increases with time?
Because that would satisfy all these criteria.
For your program this would mean introducing an extra variable Y to keep track of the area, and an extra variable t to keep track of time.
Rewriting the Runge step function, so it steps through x instead of through t.
Furthermore calculate successive values of (y, Y, t) with increasing x, instead of (x, y) with increasing t.

I can do that for you, but before I do, I'd like to know if this is useful for you, and if this is what you want.
Is it?

 

[Itosu]

If I understand correctly, what you're actually looking for is the function Y(x) which I define to be the area under the curve y(x) between x=x0 and some value of x.
Can you confirm this?

Yes.

Can I assume that x strictly increases with time?

Sample graph is like the one in the attachment. So, as you see x increases, however after some time it dramatically increases, because what I am (We are) creating is glucose-insuline model, and when glucose level increases more insuline is produced and it makes glucose level decrease.

AUC is very useful for me because it helps to model the effect of glycemic index, so if you had a little time I would be very grateful for the calculations. I need more time to think over what you just wrote :)

 

[I like Serena]

What do the axes in your graph represent?
Is it an x-t-diagram or what?

And what do your variables x, y, and Y represent?
I presume they are related to glucose level and/or insuline level?

I presume "Czas" is time in hours?
Is "Glikemia" glucose level?

 

[I like Serena]

All right.
Here's the code for the system:

dY/dx = y
dy/dx = fdy(x, y) / fdx(x, y, t)
dt/dx = 1 / fdx(x, y, t)​

        private const double H = 0.01; // This will now be the x step instead of the t step

        public double t = 0;
        public double x = 0.81;
        public double y = 0.055;
        public double Y = 0;
        // ;
        public ArrayList tList;
        public ArrayList xList;
        public ArrayList yList;
        public ArrayList YList;

        public void Start()
        {
            tList = new ArrayList();
            xList = new ArrayList();
            yList = new ArrayList();
            YList = new ArrayList();

            for (double krok = 0; krok <= 500; krok++)
            {
                xList.Add(x);
                yList.Add(y);
                YList.Add(Y);
                tList.Add(t);
                rungeStep(H, x, y, Y, t);
            }
        }

        public void rungeStep(double H, double x, double y, double Y, double t)
        {
            double K1Y = H * y;
            double K1y = H * dy(x, y) / fdx(x, y, t)
            double K1t = H * 1 / dx(x, y, t);

            double K2Y = H * (y + 0.5 * K1y);
            double K2y = H * dy(x + 0.5 * H, y + 0.5 * K1y) / dx(x + 0.5 * H, y + 0.5 * K1y, t + 0.5 * K1t);
            double K2t = H * 1 / dx(x + 0.5 * H, y + 0.5 * K1y, t + 0.5 * K1t);

            double K3Y = H * (y + 0.5 * K2y);
            double K3y = H * dy(x + 0.5 * H, y + 0.5 * K2y) / dx(x + 0.5 * H, y + 0.5 * K2y, t + 0.5 * K2t);
            double K3t = H * 1 / dx(x + 0.5 * H, y + 0.5 * K2y, t + 0.5 * K2t);

            double K4Y = H * (y + K3y);
            double K4y = H * dy(x + H, y + K3y) / dx(x + H, y + K3y, t + K3t);
            double K4t = H * 1 / dx(x + H, y + K3y, t + K3t);            double stepY = (K1Y + 2 * K2Y + 2 * K3Y + K4Y) / 6;
            double stepy = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;
            double stept = (K1t + 2 * K2t + 2 * K3t + K4t) / 6;            this.x = x + H;
            this.Y = Y + stepY;
            this.y = y + stepy;
            this.t = t + stept;
        }

 

[I like Serena]

I just realized there's another differential equation (DE) that you can use.
The advantage is that you get exactly what you got before, but you also get the area under curve as a function of time (instead of x).

The DE is:

dx/dt = fdx(x, y, t)
dy/dt = fdy(x, y)
dY/dt = y fdx(x, y, t)​

That is, just one additional evaluation to keep track of the area under curve.
It will work without any extra conditions.

 

[Itosu]

Thank you for your answers!

Could you explain in more details where dY/dt = y fdx(x, y, t) came from?

And to calculate error do I understand it correctly that I have to run the same programme with H=0.01/2 and compare results?

Thank you one more time!

 

[I like Serena]

Thank you for your answers!

Could you explain in more details where dY/dt = y fdx(x, y, t) came from?

And to calculate error do I understand it correctly that I have to run the same programme with H=0.01/2 and compare results?

Thank you one more time!

The chain rule says that:

d/dt Y(x(t)) = dY/dx dx/dt​

And since Y is the integral of y with respect to x, y is the derivative of Y with respect to x, so:

dY/dx = y​

Combining the two gives:

dY/dt = y dx/dt​

Perhaps more intuitive is thinking of rectangles with an infinitesimal width dx.
When x increases to x+dx, the area under curve increases with dY = y dx, which is the surface of an infinitesimal rectangle.
This happens between time t and t+dt.
So dY/dt = y dx/dt.
And yes, running the program again with half of the step size will give you an estimation of the error.
Of course, when you have done that, this result will be the more accurate one. :)

 

[Itosu]

Thank you. I understand everything (are you a teacher? :) but one thing:

And since Y is the integral of y with respect to x,

How do we know that? Why with respect to x, not to t?

And as far as error is concerned - I start with half smaller step size, and then substract results?

Thank you again!

 

[I like Serena]

Thank you. I understand everything (are you a teacher? :) but one thing:

How do we know that? Why with respect to x, not to t?

And as far as error is concerned - I start with half smaller step size, and then substract results?

Thank you again!

Professionally I'm not a teacher but an engineer, but I do give private lessons. :)
And as you can see, I've become a homework helper on PF.
Perhaps someday I'll be awarded a science advisor badge. ;)

You asked for the area under curve of y with respect to x.
In a previous post I have defined Y to be this area under curve (with respect to x).
Isn't that what you wanted?

As for the error.
When you calculate the result twice with the 2 different step sizes, you'll find in each case a final (t, x, y, Y).
You can interpret your error in different ways.
Since you mentioned that the area under curve is carried forward, I would interpret the error as the different areas under curve at the final time t.

Btw, I'm still interested in what you are measuring.
What do your symbols x, y and Y represent?

 

[Itosu]

I am full of admiration that you find a time to support us on the forum :)

x - glucose level
y - insuline level
Y - AUC of the glucose level function
I need AUC to calculate http://en.wikipedia.org/wiki/Glycemic_index"

So, as I need AUC under x (as in previous image), on vertical axis there is x, on horizontal time, shouldn't we estimate the integral of x?

 

[I like Serena]

I am full of admiration that you find a time to support us on the forum :)

x - glucose level
y - insuline level
Y - AUC of the glucose level function
I need AUC to calculate http://en.wikipedia.org/wiki/Glycemic_index"

Thank you!

So, as I need AUC under x (as in previous image), on vertical axis there is x, on horizontal time, shouldn't we estimate the integral of x?

Yep!
So you need a different set of DE equations.

So apparently the area under curve is $\displaystyle AUC(t) = \int_0^{t} x(t) dt$, for t = 2 hours.
This means you can use the DE set:

dx/dt = fdx(x, y, t)
dy/dt = fdy(x, y)
dAUC/dt = x​

This is much easier!

 

[Itosu]

Thank you for your answer! :)

In each step, which is equal to time interval (0.01) I calculate the value of x, so as a result I have a List<Double> of 500 x values (5 hours). But how to use them to get AUC? :)

 

[I like Serena]

Thank you for your answer! :)

In each step, which is equal to time interval (0.01) I calculate the value of x, so as a result I have a List<Double> of 500 x values (5 hours). But how to use them to get AUC? :)

I don't get it. Why would you calculate this list?

When you're solving the DE set with RK, your lists are filled with x values, y values, and AUC values.
And apparently you don't need all these lists with values!
You're only interested in the AUC at t=2 hours, aren't you?
But you do get these lists with values for free, so you can make nice plots.

 

[Itosu]

Here is my code:

 public void rungeDX(double x, double y, double t)
        {
            double K1x = H * dx(x, y, t);
            double K1y = H * dy(x, y);

            double K2x = H * dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
            double K2y = H * dy((x + 0.5 * K1x), (y + 0.5 * K1y));

            double K3x = H * dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
            double K3y = H * dy((x + 0.5 * K2x), (y + 0.5 * K2y));

            double K4x = H * dx((x + K3x), (y + K3y), t + H);
            double K4y = H * dy((x + K3x), (y + K3y));


            double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
            double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;

            this.x = x + restX;
            this.y = y + restY;
        }

and then:

            double t = 0;
            this.x = 0.81;
            this.y = 0.055;
            Solution solution = new Solution();
            for (double krok = 0; krok <= 500; krok++)
            {
                t = krok / 100;
                solution.xList.Add(x);
                solution.yList.Add(y);
                solution.time.Add(t);
                rungeDX(x, y, t);
                
            }

I need to calculate the AUC of x :)(x-glucose, y-insuline).

 

[I like Serena]

I need to calculate the AUC of x :)(x-glucose, y-insuline).

Yes. So you need to do something like:

 public void rungeDX(double x, double y, double auc, double t)
        {
            double K1x = H * dx(x, y, t);
            double K1y = H * dy(x, y);
            double K1auc = H * x;

            double K2x = H * dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
            double K2y = H * dy((x + 0.5 * K1x), (y + 0.5 * K1y));
            double K2auc = H * (x + 0.5 * K1x);

            double K3x = H * dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
            double K3y = H * dy((x + 0.5 * K2x), (y + 0.5 * K2y));
            double K3auc = H * (x + 0.5 * K2x);

            double K4x = H * dx((x + K3x), (y + K3y), t + H);
            double K4y = H * dy((x + K3x), (y + K3y));
            double K4auc = H * (x + K3x);            double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
            double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;
            double restAUC = (K1auc + 2 * K2auc + 2 * K3auc + K4auc) / 6;

            this.x = x + restX;
            this.y = y + restY;
            this.auc = auc + restAUC;
        }

and you need to declare this.auc and set it to zero before starting.

 

[Itosu]

Thank you :)
Everything is working perfectly! But if you had a moment could you explain to me the idea? We use RK to integrate x, and isn't RK designated for differential equations?

 

[I like Serena]

Thank you :)
Everything is working perfectly! But if you had a moment could you explain to me the idea? We use RK to integrate x, and isn't RK designated for differential equations?

Solving a differential equation is a more advanced form of integration.

Consider for instance the DE: dy - f(x)dx = 0
This is the same as dy/dx = f(x).
And that is the same as calculation the integral of f(x)dx.

In other words, any integral can also be written as a differential equation, where the result of the integration is the solution of the differential equation.

RK can calculate the solution for this numerically.

 

[Itosu]

Thank you!
And regarding errors - I should run the programm with step = H/2 and compare results, however on what basis? - I don't know the exact result that should be. Should I do it one more time and check if the results vary much?

 

[Itosu]

Ok, I understand that if one set of results differs much from another set (with different step) then the step should be decreased until differences won' tbe large.

One more question :)
Here we are using fixed-step. Is it difficult to implement variable step size?

Thank you in advance!

 

[I like Serena]

Thank you!
And regarding errors - I should run the programm with step = H/2 and compare results, however on what basis? - I don't know the exact result that should be. Should I do it one more time and check if the results vary much?

When you use the step H/2, the result should be about 32 times as accurate (this is what O(h⁵) means).
You can interpret this as that you have at least one extra digit that is accurate in your result.
So all the digits that are the same in both results will be presumably accurate.

Ok, I understand that if one set of results differs much from another set (with different step) then the step should be decreased until differences won' tbe large.

One more question :)
Here we are using fixed-step. Is it difficult to implement variable step size?

Thank you in advance!

No, it's not particularly difficult to implement a variable step size.

You will have to adjust your rungeDx() method to not simply update this.x, this.y, and this.auc, but these values need to be returned from the method.
That way, you can calculate each step more than once without actually executing it.

In each step you need to repeatedly calculate the result, halving the step size, until your result does not change much anymore, that is, changes less than some threshold that you defined beforehand.
With each step you can also try to repeatedly double the step size until it changes too much.

Basically, that's it.