Math Is Fun Forum

Itosu

Member

Registered: 2011-07-20

Posts: 11

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

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 am sorry, the image didn't upload. These are the equations:

Last edited by Itosu (2011-07-20 23:19:23)

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi Itosu;

Welcome to the forum.

It is not clear what type of help you require. Debugging a program is a personal activity.

Let's say your program is written correctly. Here is what I can do.
Are the DE's stiff? Are they numerically unstable? It is easy to eat up significant digits in
an iterative process such as Runge Kutta? What are the DE's ?  I can help you to solve them, maybe... Simulink could be wrong!

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

Thank you bobbym for your answer.
I uploaded the equations (sorry, I was working late at night...).

Do you have any idea how to solve them properly?

Thank you in advance!

Last edited by Itosu (2011-07-20 23:27:35)

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi;

What is dla?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

"for"

and Cg = 15000 and Ci = 15000 (divisions in dx and dy).

Last edited by Itosu (2011-07-20 23:44:47)

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

I changed solution a little, however still not what expected...

 public Solution Start()
        {
            xList = new List<Double>();
            yList = new List<Double>();
            double t = 0;

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


            }
            return solution;
        }

        public class Solution
        {
            public Solution()
            {
                time = new List<Double>();
                xList = new List<Double>();
                yList = new List<Double>();
            }
            public List<Double> time;
            public List<Double> yList;
            public List<Double> xList;
        }

        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;
        }

        public void rungeDX(double x, 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 + 0.5 * H);
            double K4y = dy((x + K3x), (y + K3y));

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

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

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi;

What do you think is the correct answer?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

I think the correct answer is what Simulink produces, as an author of the model prepared a simulation in Simulink. I can visually compare my C# graphs and those in Simulink and they are different. Could you take a look at my rungeDX method and validate?

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi;

I am far away from verifying anything right now. What I want to do is solve the DE's analytically. If I can not do that then I would like to solve them numerically. Then it will be easy to graph them. Could you telll me what U(t) and phi is? Since you are solving numerically I know you have a value for each of them or they are a function.

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

Consts are as here: (FI = PHI)

 private const int ALFA = 7600;
        private const int MU = 7200;
        private int BETA = 1430;
        private const int LAMBDA = 2470;
        private int V = 139000;
        private const double FI = 0.51; //PHI
        private const double TETA = 2.5;
        private const int QL = 8400;
        private const double H = 0.01;

and U(t) is simple function:

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

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Okay, it is a piecewise function. I did not see that.

Is Teta equal to θ ?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

Yes.

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

What I just noticed now is that both DE's are piecewise functions also.

Any idea what the interval of t is? 0 to what?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

0 to 5, with step of 0.01

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi;

What about initial conditions?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

Initial conditions:
        public double x = 0.81;
        public double y = 0.055;

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Yes, sorry I saw that and could not get to post that I did.

What kind of graph did you use, parametric , 3D, 2D?

Can you upload a picture of the graphs?

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

I think I just solved the problem:
[CODE]

        public void runge(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;

        }
[CODE]

However, now I must calculate some errors, because graphs are too little. Could you give me some piece of advice what errors and how to calculate? Thank you

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

Hi;

I am all done here. I wished I could have seen what your graph looks like so I
could compare it to mine.

You could print out a table of values and I could check them once I know mine is right.

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.

Itosu

Member

Registered: 2011-07-20

Posts: 11

Re: fourth order Runge–Kutta in C# - two differential equations

Should I paste them here?
So no errors calculating?

bobbym

bumpkin

From: Bumpkinland

Registered: 2009-04-12

Posts: 109,606

Re: fourth order Runge–Kutta in C# - two differential equations

I do not know about the errors yet. How big are the graphs in kb's.

---

In mathematics, you don't understand things. You just get used to them.
If it ain't broke, fix it until it is.
Always satisfy the Prime Directive of getting the right answer above all else.