16. august 2009 - 05:06
#5
Noget mere avanceret kode:
using System;
namespace E
{
public class Calculus
{
public delegate double Func(double x);
public static double RectangularIntegrate(Func f, double start, double finish, int nintval)
{
double eps = (finish - start) / nintval;
double sum = 0;
for(int i = 0; i < nintval; i++)
{
sum += f(start + eps * (0.5 + i));
}
return eps * sum;
}
public static double TrapezIntegrate(Func f, double start, double finish, int nintval)
{
double eps = (finish - start) / nintval;
double sum = 0.5 * f(start);
for(int i = 1; i < nintval; i++)
{
sum += f(start + eps * i);
}
sum += 0.5 * f(finish);
return eps * sum;
}
public static double SimpsonIntegrate(Func f, double start, double finish, int nintval)
{
double eps = (finish - start) / (2 * nintval);
double sum = f(start);
for(int i = 1; i <= nintval; i++)
{
sum += 4 * f(start + eps * (2 * i - 1));
}
for(int i = 1; i < nintval; i++)
{
sum += 2 * f(start + eps * 2 * i);
}
sum += f(finish);
return eps * sum / 3;
}
}
public class Program
{
private const string FMT = "{0,-30} : {1:f12} {2}";
public static double StandardNormal(double x)
{
return Math.Exp(-x*x/2)/Math.Sqrt(2*Math.PI);
}
private static double[] poly1 = { 0.000124818087,
-0.001075204047,
0.005198775019,
-0.019198292004,
0.059054035642,
-0.151968751364,
0.319152932604,
-0.531923007300,
0.797884560593 };
private static double[] poly2 = { -0.000045255659,
0.000152529290,
-0.000019538132,
-0.000676904986,
0.001390604284,
-0.000794620820,
-0.002034254874,
0.006549791214,
-0.010557625006,
0.011630447319,
-0.009279453341,
0.005353579108,
-0.002141268741,
0.000535310849,
0.999936657524 };
public static double StandardNormalAccAlgo1(double x)
{
double tmp;
double absx = Math.Abs(x);
if(absx == 0) {
tmp = 0;
} else if(absx > 6) {
tmp = 1;
} else if(absx < 2) {
double v = absx*absx/4;
tmp = 0;
for(int i = 0; i < poly1.Length; i++) tmp = tmp * v + poly1[i];
tmp = tmp * absx;
} else {
double v = absx / 2 - 2;
tmp = 0;
for(int i = 0; i < poly2.Length; i++) tmp = tmp * v + poly2[i];
}
if(x > 0)
return (tmp+1)/2;
else
return (1-tmp)/2;
}
private const double b1 = 0.319381530;
private const double b2 = -0.356563782;
private const double b3 = 1.781477937;
private const double b4 = -1.821255978;
private const double b5 = 1.330274429;
private const double p = 0.2316419;
private const double c = 0.39894228;
public static double StandardNormalAccAlgo2(double x)
{
if(x >= 0.0)
{
double t = 1.0 / ( 1.0 + p * x );
return (1.0 - c * Math.Exp( -x * x / 2.0 ) * t * ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
}
else
{
double t = 1.0 / ( 1.0 - p * x );
return ( c * Math.Exp( -x * x / 2.0 ) * t * ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
}
}
private static string Elapsed(DateTime dt1, DateTime dt2, int rep)
{
return string.Format("({0:F8} seconds)", (dt2-dt1).TotalSeconds/rep);
}
private const int REP = 100;
private static void TestRectangular(int n)
{
double res = 0;
DateTime dt1 = DateTime.Now;
for(int i = 0; i < REP; i++)
{
res = Calculus.RectangularIntegrate(StandardNormal, -1, 1, n);
}
DateTime dt2 = DateTime.Now;
Console.WriteLine(FMT, "Rectangular " + n + " slices", res, Elapsed(dt1, dt2, REP));
}
private static void TestTrapez(int n)
{
double res = 0;
DateTime dt1 = DateTime.Now;
for(int i = 0; i < REP; i++)
{
res = Calculus.TrapezIntegrate(StandardNormal, -1, 1, n);
}
DateTime dt2 = DateTime.Now;
Console.WriteLine(FMT, "Trapez " + n + " slices", res, Elapsed(dt1, dt2, REP));
}
private static void TestSimpson(int n)
{
double res = 0;
DateTime dt1 = DateTime.Now;
for(int i = 0; i < REP; i++)
{
res = Calculus.SimpsonIntegrate(StandardNormal, -1, 1, n);
}
DateTime dt2 = DateTime.Now;
Console.WriteLine(FMT, "Simpson " + n + " slices", res, Elapsed(dt1, dt2, REP));
}
private static void TestAlgo1()
{
double res = 0;
DateTime dt1 = DateTime.Now;
for(int i = 0; i < REP; i++)
{
res = StandardNormalAccAlgo1(1) - StandardNormalAccAlgo1(-1);
}
DateTime dt2 = DateTime.Now;
Console.WriteLine(FMT, "Algorithm 1", res, Elapsed(dt1, dt2, REP));
}
private static void TestAlgo2()
{
double res = 0;
DateTime dt1 = DateTime.Now;
for(int i = 0; i < REP; i++)
{
res = StandardNormalAccAlgo2(1) - StandardNormalAccAlgo1(-1);
}
DateTime dt2 = DateTime.Now;
Console.WriteLine(FMT, "Algorithm 2", res, Elapsed(dt1, dt2, REP));
}
public static void Main(string[] args)
{
Console.WriteLine("Test");
Console.WriteLine("----");
Console.WriteLine(FMT, "Correct", 0.682689492137, "");
TestRectangular(10);
TestRectangular(100);
TestRectangular(1000);
TestRectangular(10000);
TestRectangular(100000);
TestRectangular(1000000);
TestTrapez(10);
TestTrapez(100);
TestTrapez(1000);
TestTrapez(10000);
TestTrapez(100000);
TestTrapez(1000000);
TestSimpson(10);
TestSimpson(100);
TestSimpson(1000);
TestSimpson(10000);
TestSimpson(100000);
TestSimpson(1000000);
TestAlgo1();
TestAlgo2();
Console.WriteLine("Result");
Console.WriteLine("------");
Console.WriteLine(FMT, "Rectangular 1000 slices", Calculus.RectangularIntegrate(StandardNormal, -10, 1.345067, 1000), "");
Console.WriteLine(FMT, "Trapez 1000 slices", Calculus.TrapezIntegrate(StandardNormal, -10, 1.345067, 1000), "");
Console.WriteLine(FMT, "Simpson 1000 slices", Calculus.SimpsonIntegrate(StandardNormal, -10, 1.345067, 1000), "");
Console.WriteLine(FMT, "Rectangulat 1000000 slices", Calculus.RectangularIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
Console.WriteLine(FMT, "Trapez 1000000 slices", Calculus.TrapezIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
Console.WriteLine(FMT, "Simpson 1000000 slices", Calculus.SimpsonIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
Console.WriteLine(FMT, "Algorithm 1", StandardNormalAccAlgo1(1.345067), "");
Console.WriteLine(FMT, "Algorithm 2", StandardNormalAccAlgo2(1.345067), "");
Console.ReadKey();
}
}
}
