__Environment:__ [C#]

In this lesson I will show how to numerically solve algebraic and ordinary differential equations, and perform numerical integration with Simpson method. I will start with the solution of algebraic equations. The secant method is one of the simplest methods for solving algebraic equations. It is usually used as a part of a larger algorithm to improve convergence. As in any numerical algorithm, we need to check that the method is converging to a given precision in a certain number of steps. This is a precaution to avoid an infinite loop.

Our second example is a Simpson integration algorithm. The Simpson algorithm is more precise the naive integration algorithm I have used there. The basic idea of the Simpson algorithm is to sample the integrand in a number of points to get a better estimate of its variations in a given interval.

Finally, let me show a simple code for solving first order ordinary differential equations. The code uses a Runge-Kutta method. The simplest method to solve ODE is to do a Taylor expansion, which is called Euler’s method. Euler’s method approximates the solution with the series of consecutive secants. The error in Euler’s method is O(h) on every step of size h. The Runge-Kutta method has an error O(h^4) Runge-Kutta methods with a variable step size are often used in practice since they converge faster than fixed size methods.

//secant methodusing System;

class Secant

{

//declare a delegate that takes double and returns double

public delegate double Function(double x);public static void secant( int step_number,

double point1,

double point2,

Function f)

{

double p2,p1,p0,prec=.0001f; //set precision to .0001

int i;p0=f(point1);

p1=f(point2);

p2=p1-f(p1)*(p1-p0)/(f(p1)-f(p0)); //secant formula//iterate till precision goal is not met or the

// maximum number of steps is reached

for(i=0;System.Math.Abs(p2)>prec &&i<step_number;i++)

{

p0=p1;

p1=p2;

p2=p1-f(p1)*(p1-p0)/(f(p1)-f(p0));

}if(i<step_number)

Console.WriteLine(p2); //method converges

else //method does not converge

Console.WriteLine(“{0}.The method did not converge”,p2);

}

}class Demo

{ //equation f1(x)==0;public static double f1( double x)

{

return x*x*x-2*x-5;

}public static void Main()

{

Secant.secant(5,0,1,new Secant.Function(f1));

}

}//Simpson integration algorithm

using System;

//calculate the integral of f(x) between x=a and x=b

// by spliting the interval in step_number steps

class Integral

{

//declare a delegate that takes and returns double

public delegate double Function(double x);

public static double integral( Function f,

double a,

double b,

int step_number)

{

double sum=0;

double step_size=(b-a)/step_number;//Simpson algorithm samples the integrand in several

//point which significantly improves precision.

for(inti=0;i<step_number;i=i+2)

// divide the area under f(x) into step_number

// rectangles and sum their areas

sum = sum + (f(a+i*step_size)+4*f(a+(i+1)*step_size)

+ f(a+(i+2)*step_size)) *step_size/3;

return sum;

}

}class Test

{

//simple functions to be integrated

public static double f1( double x)

{

return x*x;

}public static double f2(double x)

{

return x*x*x;

}public static void Main()

{ //output the value of the integral.

Console.WriteLine(

Integral.integral(new Integral.Function(f1),

1,10,20));

}

}using System;

//fourth order Runge Kutte method for y’=f(t,y);//solve first order ode in the interval (a,b) with a given

//initial condition at x=a and fixed step h.class Runge

{

//declare a delegate that takes a double and returns

public delegate double Function(double t,double y);//double

public static void runge( double a,

double b,

double value,

double step,

Function f)

{

double t,w,k1,k2,k3,k4;

t=a;

w=value;

for(int i=0;i<(b-a)/step;i++)

{

k1=step*f(t,w);

k2=step*f(t+step/2,w+k1/2);

k3=step*f(t+step/2,w+k2/2);

k4=step*f(t+step,w+k3);

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

t=a+i*step;

Console.WriteLine(“{0} {1} “,t,w);

}

}

}class Test

{

public static double f1(double t, double y)

{

return -y+t+1;

}public static void Main()

{

Runge.runge(0,1,1,.1f,new Runge.Function(Test.f1));

}

}