Numerical Solution to Ordinary Differential Equations

Introduction



This type of equations involves the function of the independent variable and its derivatives. Example of ODE can be  dy/dx = x + y.
    Prerequisites knowledge of ordinary differential equation and how to find their exact solutions using methods like calculating particular integral and integration factor etc.

The whole idea is concentrated on tangent to problem curve.
In this tutorial we focus on: 
  • Euler's Method
  • Modified Euler's Method
  • Runga-Kutta Method( of 4th Code)    

Runga-Kutta Method( of 4th Code)


Description


This method is one of the most widely used method and again it is also an iterative method which includes the routine of well known Euler's Method.
        Suppose dy/dx = f(x,y) , y = y0 when x = x0 to find y for x = x1 = x0 + h
we have to compute 
        k1 = hf(x0, y0)
        k2 = hf(x0 + h/2, y0 + k1/2)
        k3 = hf(x0 + h/2, y0 + k2/2)
        k= hf(x0 + h/2, y0 + k3)
    k = 1/6( k1 + 2k2 + 2k3 + k4), Now the value of y1 of y for x1 = x0 +k is given by y1 = y0 + k
similarly we have to repeat for next value of x , x2 ,x... and so on

In General for xn
        k= hf(xn-1, yn-1)
        k2 = hf(xn-1 + h/2, yn-1 + k1/2)
        k3 = hf(xn-1 + h/2, yn-1 + k2/2)
        k4 = hf(xn-1 + h/2, yn-1 + k3)
    k = 1/6( k1 + 2k2 + 2k+ k4),   yn = yn-1 + k


  • Programming implementation of RK4 in C++
#include<iostream>
using namespace std;

double f(double x,double y){
    return x*y + 1;
}
double f1(double x, double y){
    return 3*x + y/2;
}
int main(){
    //Suppose dy/dx = f(x,y) and y = y0 when x=x0 
    double x0,y0;
    double xn,yn;
    double k1,k2,k3,k4, k;
    cout<<"Provide initial (x0,y0) whose differential equ is defined:";
     cin>>x0>>y0;
     cout<<"Enter the value of x for which you need to find out the y=?:";
     cin>>xn;

     double h = (xn-x0)/5;
     yn = y0;
     for(double xi = x0 ; xi < xn ; xi += h){
         k1 = h*f(xi , yn);
         k2 = h*f(xi + h/2, yn + k1/2);
         k3 = h*f(xi + h/2, yn + k2/2);
         k4 = h*f(xi + h/2, yn + k3);

         k = (k1 + 2*k2 + 2*k3 + k4)/6;

         yn = yn+k;
     }

     printf("Value of y at x=%f is %0.4f when h: %f\n",xn , yn, h);
    return 0;
}
  • Code Execution of RK4 for dy/dx = xy +1 , y(0)=1 at x = 0.1
Provide initial (x0,y0) whose differential equ is defined:0 1
Enter the value of x for which you need to find out the y=?:0.1
Value of y at x=0.100000 is 1.1052 when h: 0.020000
    Note: check the value of h, calculate the accuracy of the method and compare with other methods of finding the solution of the ODE.


Modified Euler's Method


Description

Modified Euler's Method is just the extended version of Euler's Method to increase the accuracy of the approximate solution.
    It's recommended for you to go through Euler's Method.
The approximate value of y1 of y for x = x1 = x0 +h by Euler's method is given by 
                y1 = y+ hf(x0, y0)
so in edition to attend more accuracy we repeat the below formula, 
                y1(1) = y0 + h/2[ f(x0, y0) + f(x1, y1) ]
                y1(2) = y0 + h/2[ f(x0, y0) + f(x1, y1(1))]
                ...
in pan and paper process continues untill two value coincides now we have value of y1
                y2 = y1 + hf(x1, y1)
repeat the above process for y2
                y2(1) = y1 + h/2[ f(x1, y1) + f(x2, y2) ]
                y2(2) = y1 + h/2[ f(x1, y1) + f(x2, y2(1))]
                ...
  • Correctness

  • Programming implementation of modified Euler's Method in C++
#include<iostream>
#include<math.h>
using namespace std;

double f(double x,double y){
    return x+y;  // dy/dx = x + y
}
double f1(double x,double y){
    return log2(x+y);
}
int main(){
    // COde for solving Ordinary differential equations
    //     dy/dx = f(x,y) 
    //Define your equation

    //Given y(x0) = y0
    //Need to find y(xn) = yn
    double x0,y0;
    double xn,yn;

     cout<<"Provide initial (x0,y0) whose differential equ is defined:";
     cin>>x0>>y0;
     cout<<"Enter the value of x for which you need to find out the y=?:";
     cin>>xn;

     double h = (xn-x0)/5;
     yn = y0;
     for(double xi = x0;xi < xn; xi += h){
         cout<<xi<<" "<<yn<<endl;
         double YnPrevious = yn;

         yn = yn + h*f(xi , yn);
         for(int i=0;i<5;i++){
             yn = YnPrevious + h*( f(xi,YnPrevious) + f(xi+h,yn))/2;
         }
     //  Yn = Yn-1 + h*f(Xn-1,Yn-1) 
     }

     printf("Value of y at x=%f is %0.4f when h: %f\n",xn , yn, h);
    return 0;
}
  • Code Execution for dy/dx = x + y ,y(0) =1 ,for x=0.3
Provide initial (x0,y0) whose differential equ is defined:0 1
Enter the value of x for which you need to find out the y=?:0.3
0 1
0.06 1.06371
0.12 1.13507
0.18 1.21456
0.24 1.30268
Value of y at x=0.300000 is 1.4000 when h: 0.060000
    Note: In code, h is defined constant depending upon interval size, but you can change as h cause accuracy difference, try out the code for other functions. 


Euler's Method


Description


Consider a defferential equation dy/dx = f(x,y) with initial condition y = y0 when x = x0, suppose we want to find y for x = xn.
    First we divide the interval [x0,xn] into n equal parts. x0, x1 = x0 +h , x2 = x0 + 2h, .. , xn where h is small. Using the equation of tangent at point (x0, y0)
            y - y0 = dy/dx * (x - x0)  where dy/dx = f(x0,y0)

In the small interval [x0, x1] or [x0, x0+h]  we approximate the solution curve by the tangent, So the approximate value of y1 at point y1 can be given by
            y1 - y0 = f(x0, y0)( x0 + h - x0)
            y1 = y0 + hf(x0, y0)
Similarly
            y2 = y+ hf(x1, y1)
            ...
            yn = yn-1 + hf(xn-1, yn-1)
            


In the Euler's method we approximate the solution curve by tangent in each subinterval by a sequence of short lines and therefore unless h is very small the error is bound to be quite significant. 
  • Programming implementation of Euler's Method in C++
#include<iostream>
#include<math.h>
using namespace std;

double f(double x,double y){
    return x+y;  // dy/dx = x + y
}
double f1(double x,double y){
    return log2(x+y);
}
int main(){
    // COde for solving Ordinary differential equations
    //     dy/dx = f(x,y) 
    //Define your equation

    //Given y(x0) = y0
    //Need to find y(xn) = yn
    double x0,y0;
    double xn,yn;

     cout<<"Provide initial (x0,y0) whose differential equ is defined:";
     cin>>x0>>y0;
     cout<<"Enter the value of x for which you need to find out the y=?:";
     cin>>xn;

     double h = (xn-x0)/5;
     yn = y0;
     for(double xi = x0;xi < xn; xi += h){
         cout<<yn<<endl;
         yn = yn + h*f(xi , yn);
     //  Yn = Yn-1 + h*f(Xn-1,Yn-1) 
     }

     printf("Value of y at x=%f is %0.4f when h: %f\n",xn , yn, h);
    return 0;
}
  • Code execution for dy/dx = x + y ,y(0) = 1 ,finding y for x = 0.1
Provide initial (x0,y0) whose differential equ is defined:0 1
Enter the value of x for which you need to find out the y=?:0.1
1
1.02
1.0408
1.06242
1.08486
Value of y at x=0.100000 is 1.1082 when h: 0.020000

    Note: In code their h is defined constant, change it if needed as accuracy completely dependent on the value of h.