# Gill’s 4th Order Method to solve Differential Equations

Given the following inputs:

• An ordinary differential equation that defines the value of dy/dx in the form x and y. • Initial value of y, i.e., y(0). The task is to find the value of unknown function y at a given point x, i.e. y(x).

Examples:

Input: x0 = 0, y = 3.0, x = 5.0, h = 0.2
Output: y(x) = 3.410426

Input: x0 = 0, y = 1, x = 3, h = 0.3
Output: y(x) = 1.669395

## Recommended: Please try your approach on {IDE} first, before moving on to the solution.

Approach:
Gill’s method is used to find an approximate value of y for a given x. Below is the formula used to compute next value yn+1 from previous value yn.
Therefore:

yn+1 = value of y at (x = n + 1)
yn = value of y at (x = n)
where
0 ≤ n ≤ (x - x0)/h
h is step height
xn+1 = x0 + h


The essential formula to compute the value of y(n+1):     The formula basically computes the next value yn+1 using current yn:

• K1 is the increment based on the slope at the beginning of the interval, using y.
• K2 is the increment based on the slope, using • K3 is the increment based on the slope, using • K4 is the increment based on the slope, using The method is a fourth-order method, meaning that the local truncation error is of the order of O(h5).

Below is the implementation of the above approach:

## C++

 // C++ program to implement Gill's method     #include  using namespace std;     // A sample differential equation  // "dy/dx = (x - y)/2"  float dydx(float x, float y)  {      return (x - y) / 2;  }     // Finds value of y for a given x  // using step size h and initial  // value y0 at x0  float Gill(float x0, float y0,          float x, float h)  {      // Count number of iterations      // using step size or height h      int n = (int)((x - x0) / h);         // Value of K_i      float k1, k2, k3, k4;         // Initial value of y(0)      float y = y0;         // Iterate for number of iteration      for (int i = 1; i <= n; i++) {             // Apply Gill's Formulas to          // find next value of y             // Value of K1          k1 = h * dydx(x0, y);             // Value of K2          k2 = h * dydx(x0 + 0.5 * h,                      y + 0.5 * k1);             // Value of K3          k3 = h * dydx(x0 + 0.5 * h,                      y + 0.5 * (-1 + sqrt(2)) * k1                          + k2 * (1 - 0.5 * sqrt(2)));             // Value of K4          k4 = h * dydx(x0 + h,                      y - (0.5 * sqrt(2)) * k2                          + k3 * (1 + 0.5 * sqrt(2)));             // Find the next value of y(n+1)          // using y(n) and values of K in          // the above steps          y = y + (1.0 / 6) * (k1 + (2 - sqrt(2)) * k2                               + (2 + sqrt(2)) * k3 + k4);             // Update next value of x          x0 = x0 + h;      }         // Return the final value of dy/dx      return y;  }     // Driver Code  int main()  {      float x0 = 0, y = 3.0,          x = 5.0, h = 0.2;         printf("y(x) = %.6f",          Gill(x0, y, x, h));      return 0;  }

## Java

 // Java program to implement Gill's method  class GFG{     // A sample differential equation  // "dy/dx = (x - y)/2"  static double dydx(double x, double y)  {      return (x - y) / 2;  }     // Finds value of y for a given x  // using step size h and initial  // value y0 at x0  static double Gill(double x0, double y0,                     double x, double h)  {             // Count number of iterations      // using step size or height h      int n = (int)((x - x0) / h);         // Value of K_i      double k1, k2, k3, k4;         // Initial value of y(0)      double y = y0;         // Iterate for number of iteration      for(int i = 1; i <= n; i++)      {                    // Apply Gill's Formulas to         // find next value of y                   // Value of K1         k1 = h * dydx(x0, y);                   // Value of K2         k2 = h * dydx(x0 + 0.5 * h,                        y + 0.5 * k1);                   // Value of K3         k3 = h * dydx(x0 + 0.5 * h,                        y + 0.5 * (-1 + Math.sqrt(2)) *                       k1 + k2 * (1 - 0.5 * Math.sqrt(2)));                   // Value of K4         k4 = h * dydx(x0 + h,                        y - (0.5 * Math.sqrt(2)) *                       k2 + k3 * (1 + 0.5 * Math.sqrt(2)));                   // Find the next value of y(n+1)         // using y(n) and values of K in         // the above steps         y = y + (1.0 / 6) * (k1 + (2 - Math.sqrt(2)) *                               k2 + (2 + Math.sqrt(2)) *                              k3 + k4);                   // Update next value of x         x0 = x0 + h;      }         // Return the final value of dy/dx      return y;  }     // Driver Code  public static void main(String[] args)  {      double x0 = 0, y = 3.0,              x = 5.0, h = 0.2;         System.out.printf("y(x) = %.6f", Gill(x0, y, x, h));  }  }     // This code is contributed by Amit Katiyar

## Python3

 # Python3 program to implement Gill's method  from math import sqrt     # A sample differential equation  # "dy/dx = (x - y)/2"  def dydx(x, y):      return (x - y) / 2    # Finds value of y for a given x  # using step size h and initial  # value y0 at x0  def Gill(x0, y0, x, h):             # Count number of iterations      # using step size or height h      n = ((x - x0) / h)         # Initial value of y(0)      y = y0         # Iterate for number of iteration      for i in range(1, int(n + 1), 1):                     # Apply Gill's Formulas to          # find next value of y             # Value of K1          k1 = h * dydx(x0, y)             # Value of K2          k2 = h * dydx(x0 + 0.5 * h,                          y + 0.5 * k1)             # Value of K3          k3 = h * dydx(x0 + 0.5 * h,                         y + 0.5 * (-1 + sqrt(2)) *                        k1 + k2 * (1 - 0.5 * sqrt(2)))             # Value of K4          k4 = h * dydx(x0 + h, y - (0.5 * sqrt(2)) *                     k2 + k3 * (1 + 0.5 * sqrt(2)))             # Find the next value of y(n+1)          # using y(n) and values of K in          # the above steps          y = y + (1 / 6) * (k1 + (2 - sqrt(2)) *                             k2 + (2 + sqrt(2)) *                            k3 + k4)             # Update next value of x          x0 = x0 + h         # Return the final value of dy/dx      return y     # Driver Code  if __name__ == '__main__':             x0 = 0     y = 3.0     x = 5.0     h = 0.2        print("y(x) =", round(Gill(x0, y, x, h), 6))     # This code is contributed by Surendra_Gangwar

## C#

 // C# program to implement Gill's method  using System;     class GFG{      // A sample differential equation  // "dy/dx = (x - y)/2"  static double dydx(double x, double y)  {      return (x - y) / 2;  }      // Finds value of y for a given x  // using step size h and initial  // value y0 at x0  static double Gill(double x0, double y0,                     double x, double h)  {              // Count number of iterations      // using step size or height h      int n = (int)((x - x0) / h);          // Value of K_i      double k1, k2, k3, k4;          // Initial value of y(0)      double y = y0;          // Iterate for number of iteration      for(int i = 1; i <= n; i++)      {                     // Apply Gill's Formulas to         // find next value of y                    // Value of K1         k1 = h * dydx(x0, y);                    // Value of K2         k2 = h * dydx(x0 + 0.5 * h,                        y + 0.5 * k1);                    // Value of K3         k3 = h * dydx(x0 + 0.5 * h,                        y + 0.5 * (-1 + Math.Sqrt(2)) *                       k1 + k2 * (1 - 0.5 * Math.Sqrt(2)));                    // Value of K4         k4 = h * dydx(x0 + h,                        y - (0.5 * Math.Sqrt(2)) *                       k2 + k3 * (1 + 0.5 * Math.Sqrt(2)));                    // Find the next value of y(n+1)         // using y(n) and values of K in         // the above steps         y = y + (1.0 / 6) * (k1 + (2 - Math.Sqrt(2)) *                               k2 + (2 + Math.Sqrt(2)) *                              k3 + k4);                    // Update next value of x         x0 = x0 + h;      }          // Return the final value of dy/dx      return y;  }      // Driver Code  public static void Main(String[] args)  {      double x0 = 0, y = 3.0,              x = 5.0, h = 0.2;          Console.Write("y(x) = {0:F6}", Gill(x0, y, x, h));  }  }     // This code is contributed by Amit Katiyar

Output:

y(x) = 3.410426


Attention reader! Don’t stop learning now. Get hold of all the important DSA concepts with the DSA Self Paced Course at a student-friendly price and become industry ready.

My Personal Notes arrow_drop_up Check out this Author's contributed articles.

If you like GeeksforGeeks and would like to contribute, you can also write an article using contribute.geeksforgeeks.org or mail your article to contribute@geeksforgeeks.org. See your article appearing on the GeeksforGeeks main page and help other Geeks.

Please Improve this article if you find anything incorrect by clicking on the "Improve Article" button below.

Article Tags :
Practice Tags :

Be the First to upvote.

Please write to us at contribute@geeksforgeeks.org to report any issue with the above content.