InterpolationPolynom.Cs – класс, реализующий вычисление интерполяционного полинома.

using System;

using System.Linq;

using NUnit.Framework;

using NumericAnalysis.FunctionHelpers;

namespace NumericAnalysis.Interpolation {

public class InterpolationPolynom {

private readonly double[] x;

private readonly double[] f;

private readonly int n;

public double Calc(double arg) {

double S = f[n];

for (int i = n - 1; i >= 0; i--) {

S *= (arg - x[i]);

S += f[i];

}

return S;

}

public InterpolationPolynom(Calc function, double[] x, double epsilon ) {

this.x = x.ToArray();

this.n = x.Length - 1;

this.f = x.Select(_x => function(_x, epsilon)).ToArray();

for (int k = 1; k <= n; k++) {

for (int i = n; i >= k; i--) {

this.f[i] = (this.f[i] - this.f[i - 1]) / (x[i] - x[i - k]);

}

}

}

public InterpolationPolynom(double[] x, double[] f) {

this.x = x;

this.n = f.Length - 1;

this.f = f.ToArray();

for(int k = 1; k<= n; k++) {

for (int i = n; i >= k; i--) {

this.f[i] = (this.f[i] - this.f[i - 1])/(x[i] - x[i - k]);

}

}

}

}

}

InverseFunctionHelper.Cs – класс, реализующий вычисление обратной функции.

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using NumericAnalysis.FunctionHelpers;

namespace NumericAnalysis.InverseFunction {

class InverseFunctionHelper {

public double calc (double F0, double x0, double epsilon) {

double z0 = x0;//, z1 = x0;

double z2;

Calc _g = BesselFunctionHelper.Instance.Calc;

Func<double, double> g = d => _g(d, epsilon) - F0;

Calc _g_ = BesselFunctionDerivativeHelper.Instance.Calc;

Func<double, double> g_ = d => _g_(d, epsilon);

//double gz1 = g(z1);

//double g_z1 = g_(z1);

double gz0 = g(z0);

double g_z0 = g_(z0);

if (Math.Abs(g_z0 - 0.0) < 0.000001)

return z0;

double z1 = z0 - (gz0/g_z0);

do {

double gz1 = g(z1);

double g_z1 = g_(z1);

gz0 = g(z0);

g_z0 = g_(z0);

z2 = z1 - (gz1 / g_z1) + (gz1*gz1/(gz1 - gz0))*(2/g_z1 + 1/g_z0 - 3*(z1 - z0)/(gz1 - gz0));

if(Math.Abs(z2 - z1) < epsilon) break;

z0 = z1;

z1 = z2;

} while (true);

return z2;

}

}

}

Program.cs – проводит вычисления, необходимые для выполнения задания.

using System;

using System.Collections.Generic;

using System.Diagnostics;

using System.Linq;

using System.Text;

using NumericAnalysis.FunctionHelpers;

using NumericAnalysis.Integral;

using NumericAnalysis.Interpolation;

using NumericAnalysis.InverseFunction;

namespace NumericAnalysis {

class Program {

static void Main(string[] args)

{

Main2();

}

static void Main2() {

var helper = new InverseFunctionHelper();

//double xi = 0.15;

//double Fi = BesselFunctionHelper.Instance.Calc(xi, 0.000001);

int N = 11;

double[] xi = new double[N];

double a = 0.0;

double h = 0.3;

for (int i = 0; i < 11; i++ ) {

xi[i] = a + i*h;

}

double [] fi = new double[N];

double f0 = BesselFunctionHelper.Instance.Calc(xi[0], 0.000001);

double fN = BesselFunctionHelper.Instance.Calc(xi[N - 1], 0.000001);

for (int i = 0; i < 11; i++ ) {

fi[i] = f0 + i*(fN - f0)/N;

}

for (int i = 0; i < N; i++) {

var zi = helper.calc(fi[i], xi[i], 0.000001);

Console.WriteLine("F{0} = {1}, x{0} = {2}, J{0} = {3}", i, fi[i], zi,

BesselFunctionHelper.Instance.Calc(zi, 0.000001));

}

Console.ReadLine();

}

static void Main1(string[] args) {

Stopwatch sw = Stopwatch.StartNew();

//BesselFunctionHelper besselFunctionHelper = new BesselFunctionHelper();

const double a = 0;

const double b = 3;

const double step = 0.3;

const double epsilon = 0.000001;

//

double[] x = Enumerable.Range(0, Convert.ToInt32(((b - a)/step)) + 1).Select(i => a + i*step).ToArray();

double[] _f = BesselFunctionHelper.Instance.CalcMany(x, epsilon);

for(int i = 0; i < x.Length; i++)

{

Console.WriteLine("{0:F2} {1:F6}", x[i], _f[i]);

}

//

const int n = 40;

Random r= new Random();

double[] x2 = Enumerable.Range(0, n).Select(i => a + r.NextDouble()*(b - a)).OrderBy(_x2 => _x2).ToArray();

//

InterpolationPolynom polynom = new InterpolationPolynom(BesselFunctionHelper.Instance.Calc, x, epsilon);//new InterpolationPolynom(x, f);

double[] f2 = x2.Select(polynom.Calc).ToArray();

double[] f2_2 = x2.Select(_x2 => BesselFunctionHelper.Instance.Calc(_x2, epsilon)).ToArray();// besselFunctionHelper.CalcMany(x2, epsilon);

double[] _epsilon = new double[f2.Length];

for (int i = 0; i < f2.Length; i++ ) {

_epsilon[i] = Math.Abs(f2_2[i] - f2[i]);

}

double max_epsilon = _epsilon.Max();

for(int i = 0; i <x.Length; i++) {

Console.WriteLine("{0} {1} {2} {3}",

String.Format("{0:F2}", x2[i]),

String.Format("{1}{0:F6}", f2[i], f2[i] > 0 ? " " : ""),

String.Format("{1}{0:F6}", f2_2[i], f2_2[i] > 0 ? " " : ""),

String.Format("{0:F6}", _epsilon[i])

);

}

Console.WriteLine();

Console.WriteLine("Max epsilon: ");

Console.WriteLine(max_epsilon);

//

Console.WriteLine("Trapezoidal method");

foreach(double _x in x)

{

int N;

double f = BesselFunctionHelper.Instance.CalcWithIntegral(SumAlgorithm.Trapezoidal, _x, 0.0001, out N);

Console.WriteLine("J_o({0:f2}) = {1:F4}; N = {2}", _x, f, N);

}

//

Console.WriteLine("Simpson method");

foreach (double _x in x) {

int N;

double f = BesselFunctionHelper.Instance.CalcWithIntegral(SumAlgorithm.Simpson, _x, 0.0001, out N);

Console.WriteLine("J_o({0:f2}) = {1:F4}; N = {2}", _x, f, N);

}

//

Console.WriteLine("Gauss method");

foreach (double _x in x) {

int N;

double f = BesselFunctionHelper.Instance.CalcWithIntegral(SumAlgorithm.Gauss, _x, 0.0001, out N);

Console.WriteLine("J_o({0:f2}) = {1:F4}; N = {2}", _x, f, N);

}

//

Console.WriteLine("Rectangle method");

foreach (double _x in x) {

int N;

double f = BesselFunctionHelper.Instance.CalcWithIntegral(SumAlgorithm.Rectangle, _x, 0.0001, out N);

Console.WriteLine("J_o({0:f2}) = {1:F4}; N = {2}", _x, f, N);

}

//

sw.Stop();

Console.WriteLine("Time elapsed: {0} ms.", sw.ElapsedMilliseconds);

}

}

}