using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
namespace NumMethods2
{
public partial class Form1 : Form
{
int seriesSize = 0;
double[] dataX = new double[1000000];
double[] dataY = new double[1000000];
double ans1, ans2, ans3;
class Mtr
{
double[,] data;
int n, k;
public double this[int i, int j]
{
get
{
return data[i, j];
}
set
{
data[i, j] = value;
}
}
public Mtr(int n, int k)
{
this.n = n;
this.k = k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = 0;
}
public Mtr(Mtr c)
{
this.n = c.n;
this.k = c.k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = c.data[i, j];
}
public int getColumns()
{
return k;
}
public int getRows()
{
return n;
}
public void makeIdentity()
{
if (n != k)
{
Console.WriteLine("MAKE IDENTITY ERROR: N != K");
return;
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
{
data[i, j] = ((i == j) ? 1 : 0);
}
}
public Mtr(double[,] d)
{
n = d.GetLength(0);
k = d.GetLength(1);
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = d[i, j];
}
public void setValue(double v, int i, int j)
{
data[i, j] = v;
}
public double getValue(int i, int j)
{
return data[i, j];
}
static public Mtr operator *(Mtr inp1, Mtr inp2)
{
if (inp1.k != inp2.n)
{
Console.WriteLine("MATRIX MULT ERROR");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
double tmp = 0;
for (int c = 0; c < inp1.k; ++c)
tmp += inp1.data[i, c] * inp2.data[c, j];
rs.setValue(tmp, i, j);
}
}
return rs;
}
static public Mtr operator *(Mtr inp1, double alpha)
{
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] * alpha, i, j);
}
}
return rs;
}
static public Mtr operator -(Mtr inp1, double alpha)
{
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - alpha, i, j);
}
}
return rs;
}
static public Mtr operator -(Mtr inp1, Mtr inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator- ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - inp2.data[i, j], i, j);
}
}
return rs;
}
static public Mtr operator +(Mtr inp1, Mtr inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator+ ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] + inp2.data[i, j], i, j);
}
}
return rs;
}
/// Tensor multiplication of vectors
static public Mtr operator ^(Mtr inp1, Mtr inp2)
{
if (inp1.k != inp2.n || inp1.n != inp2.k)
{
Console.WriteLine("Operator^ ERROR: Matr size not equal");
return new Mtr(1, 1);
}
Mtr rs = new Mtr(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
rs.setValue(inp1.data[i, 0] * inp2.data[0, j], i, j);
}
}
return rs;
}
public Mtr getTranspozed()
{
Mtr ret = new Mtr(k, n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret.setValue(data[i, j], j, i);
}
}
return ret;
}
public void transpoze()
{
double[,] ret = new double[k, n];
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret[j, i] = data[i, j];
}
}
data = ret;
}
public Mtr getSubMatrix(int si, int sj, int ei, int ej)
{
if (si >= n || sj >= k || ei >= n || ej >= k ||
si < 0 || sj < 0 || ei < 0 || ej < 0)
{
Console.WriteLine("getSubMatrix RANGE ERROR");
return new Mtr(1, 1);
}
Mtr retMtr = new Mtr(ei - si + 1, ej - sj + 1);
int fi = 0;
int fj = 0;
for (int i = si; i <= ei; ++i)
{
fj = 0;
for (int j = sj; j <= ej; ++j)
{
retMtr.setValue(data[i, j], fi, fj);
++fj;
}
++fi;
}
return retMtr;
}
public bool isZeroMatrix()
{
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
if (data[i, j] != 0)
return false;
return true;
}
public void printMatrix()
{
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
Console.Write("{0, 3:F4} ", data[i, j]);
}
Console.WriteLine();
}
Console.WriteLine(); Console.WriteLine(); Console.WriteLine();
}
public double getMaximumAboveDiag(ref int ii, ref int jj)
{
double val = -1000 * 1000 * 1000;
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < k; ++j)
{
if (val < Math.Abs(data[i, j]))
{
val = Math.Abs(data[i, j]);
ii = i;
jj = j;
}
}
}
return data[ii, jj];
}
}
Mtr func(double x, double y, double z)
{
Mtr ans = new Mtr(3, 1);
ans[0, 0] = 2 * x * x - y - 1;
ans[1, 0] = 2 * y * y - 3 * z - 2;
ans[2, 0] = 2 * z * z - x - 1;
return ans;
}
const double delta = 1e-9;
Mtr getJacobianAt(double x, double y, double z)
{
Mtr fnc = func(x, y, z);
Mtr dx = (func(x + delta, y, z) - fnc) * (1.0 / delta);
Mtr dy = (func(x, y + delta, z) - fnc) * (1.0 / delta);
Mtr dz = (func(x, y, z + delta) - fnc) * (1.0 / delta);
Mtr jacobian = new Mtr(3, 3);
jacobian[0, 0] = dx[0, 0];
jacobian[1, 0] = dx[1, 0];
jacobian[2, 0] = dx[2, 0];
jacobian[0, 1] = dy[0, 0];
jacobian[1, 1] = dy[1, 0];
jacobian[2, 1] = dy[2, 0];
jacobian[0, 2] = dz[0, 0];
jacobian[1, 2] = dz[1, 0];
jacobian[2, 2] = dz[2, 0];
return jacobian;
}
void LuDecomposition(Mtr A, ref Mtr L, ref Mtr U)
{
if (A.getRows() != A.getColumns())
return;
int n = A.getRows();
L.makeIdentity();
for (int i = 0; i < n; ++i)
{
for (int j = i; j < n; ++j)
{
double temp = A.getValue(i, j);
for (int k = 0; k < i; ++k)
{
temp = temp - U.getValue(k, j) * L.getValue(i, k);
}
U.setValue(temp, i, j);
}
for (int j = i + 1; j < n; ++j)
{
double temp = A.getValue(j, i);
for (int k = 0; k < i; ++k)
{
temp = temp - U.getValue(k, i) * L.getValue(j, k);
}
temp = temp / U.getValue(i, i);
L.setValue(temp, j, i);
}
}
}
Mtr inversedLuSolver(Mtr A, Mtr b)
{
int n = A.getColumns();
Mtr L = new Mtr(n, n);
Mtr U = new Mtr(n, n);
LuDecomposition(A, ref L, ref U);
Mtr inversed = new Mtr(n, n);
for (int j = n - 1; j >= 0; --j)
{
double tmp = 1;
for (int k = j + 1; k < n; ++k)
tmp -= U.getValue(j, k) * inversed.getValue(k, j);
tmp /= U.getValue(j, j);
inversed.setValue(tmp, j, j);
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += U.getValue(i, k) * inversed.getValue(k, j);
tmp = -tmp / U.getValue(i, i);
inversed.setValue(tmp, i, j);
}
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += inversed.getValue(j, k) * L.getValue(k, i);
inversed.setValue(-tmp, j, i);
}
}
Mtr answer = inversed * b;
return answer;
}
Mtr getInversed(Mtr A)
{
int n = A.getColumns();
Mtr L = new Mtr(n, n);
Mtr U = new Mtr(n, n);
LuDecomposition(A, ref L, ref U);
Mtr inversed = new Mtr(n, n);
for (int j = n - 1; j >= 0; --j)
{
double tmp = 1;
for (int k = j + 1; k < n; ++k)
tmp -= U.getValue(j, k) * inversed.getValue(k, j);
tmp /= U.getValue(j, j);
inversed.setValue(tmp, j, j);
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += U.getValue(i, k) * inversed.getValue(k, j);
tmp = -tmp / U.getValue(i, i);
inversed.setValue(tmp, i, j);
}
for (int i = j - 1; i >= 0; --i)
{
tmp = 0;
for (int k = i + 1; k < n; ++k)
tmp += inversed.getValue(j, k) * L.getValue(k, i);
inversed.setValue(-tmp, j, i);
}
}
return inversed;
}
// Метод ньютона
// Суть метода: есть начальное приближение, пока не выполнено условие остановки
// разность между предыдущим значением и полученным сейчас меньше эпсилон
// f(Xn+1) < e, вычисляем новое приближение
// Критерий сходимости:
void getNewtone()
{
Mtr answer = new Mtr(3, 1);
int fix = 2;
answer[0, 0] = -0.980616362846145;
answer[1, 0] = 0.923216902163204;
answer[2, 0] = -0.098447034373452;
if (fix == 0)
answer[0, 0] = -10;
if (fix == 1)
answer[1, 0] = 10;
if (fix == 2)
answer[2, 0] = -10;
Mtr jcb = getJacobianAt(answer[0, 0], answer[1, 0], answer[2, 0]);
jcb = getInversed(jcb);
for (int i = 0; i < 140; ++i)
{
dataY[seriesSize++] = answer[fix, 0];
answer[fix, 0] = answer[fix, 0] - (jcb * func(answer[0, 0], answer[1, 0], answer[2, 0]))[fix, 0];
jcb = getJacobianAt(answer[0, 0], answer[1, 0], answer[2, 0]);
jcb = getInversed(jcb);
}
Mtr vl = func(answer[0, 0], answer[1, 0], answer[2, 0]);
ans1 = answer[0, 0];
ans2 = answer[1, 0];
ans3 = answer[2, 0];
}
// Зейделевский метод просто итерации
// Суть метода:
// Критерий сходимости:
void getZeidelMethod()
{
Mtr answer = new Mtr(3, 1);
int fix = 2;
answer[0, 0] = -0.980616362846145;
answer[1, 0] = 0.923216902163204;
answer[2, 0] = -0.098447034373452;
if (fix == 0)
answer[0, 0] = -10;
if (fix == 1)
answer[1, 0] = 10;
if (fix == 2)
answer[2, 0] = -10;
Mtr jcb = getJacobianAt(answer[0, 0], answer[1, 0], answer[2, 0]);
jcb = getInversed(jcb);
for (int i = 0; i < 1540; ++i)
{
dataY[seriesSize++] = answer[fix, 0];
answer[fix, 0] = answer[fix, 0] - (jcb * func(answer[0, 0], answer[1, 0], answer[2, 0]))[fix, 0];
}
Mtr vl = func(answer[0, 0], answer[1, 0], answer[2, 0]);
ans1 = answer[0, 0];
ans2 = answer[1, 0];
ans3 = answer[2, 0];
}
double func2(double x)
{
return -x * x * x + x * 2 + 10 + x * x;
}
double func3(double x)
{
return x - func2(x);
}
void getChordMethod()
{
double answer;
double a = -200;
double b = 200;
double c = 0, oldc = 10000;
while (Math.Abs(c - oldc) > 1e-13)
{
oldc = c;
c = a - (func2(a) * (b - a)) / (func2(b) - func2(a));
if (func2(c) * func2(a) < 0)
{
b = c;
}
else
{
a = c;
}
dataX[seriesSize] = seriesSize;
dataY[seriesSize++] = a;
}
Console.WriteLine("Метод хорд: ");
if (Math.Abs(func2(a)) < Math.Abs(func2(b)))
ans1 = a;
else
ans1 = b;
}
void getIterMethod()
{
double x = 10, df;
df = (func2(x + delta) - func2(x - delta)) / (2 * delta);
for (int i = 0; i < 200; ++i)
{
dataX[seriesSize] = i;
dataY[seriesSize++] = x;
x = x - func2(x) / df;
}
ans1 = x;
}
void getVegseinMethod()
{
double eps = 1e-8;
double x1 = 0, x2 = 0, x0 = 1.7, x1s = 0, x2s = 0, x0s = 0;
x1 = func3(x0);
x0s = x0;
x1s = x1;
for (int i = 0; i < 20; ++i)
{
double df = Math.Abs((func3(x1s + delta) - func3(x1s - delta)) / (2 * delta));
x2 = func3(x1s);
dataX[seriesSize] = i;
dataY[seriesSize++] = x2;
if (Math.Abs(x2 - x1s) > eps * (1 - df) / df)
{
x2s = (x2 * x0s - x1 * x1s) / (x2 + x0s - x1 - x1s);
}
x0 = x1;
x0s = x1s;
x1 = x2;
x1s = x2s;
}
ans1 = x2;
}
double[] xarr = new double[] { 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4 };
double[] yarr = new double[] { 0, 0.45, 0.63, 0.77, 0.89, -1.0, 1.1, 1.18};
// double[] xarr = new double[] { -2, 1, 4, 7 };
// double[] yarr = new double[] { 0, -1, 0, -2 };
double[] factors = new double[] { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880 };
double getDelta(int m, int n, double[] yarr)
{
if (m == 1)
return yarr[n + 1] - yarr[n];
else
return getDelta(m - 1, n + 1, yarr) - getDelta(m - 1, n, yarr);
}
double fncForInversed(double x)
{
return Math.Log(x * Math.Log(1 + x)) / x ;
}
// Обратная интерполяция
// Суть метода:
// Критерий сходимости:
double findRootByInversedInterpolation(double y0)
{
double x0 = 1.0;
double h = 0.1;
int n = 5;
double[] yarr = new double[n];
for (int i = 0; i < n; ++i)
yarr[i] = fncForInversed(x0 + i * h);
double q = (y0 - yarr[0]) / getDelta(1, 0, yarr);
for (int it = 0; it < 10; ++it)
{
int factor = 1;
double nq = (y0 - yarr[0]) / getDelta(1, 0, yarr);
double mult = q;
for (int i = 2; i < n; ++i)
{
factor *= i;
mult *= (q - i + 1);
nq -= getDelta(i, 0, yarr) / (getDelta(1, 0, yarr) * factor) * mult;
}
q = nq;
}
double ans = x0 + q * h;
double funcVal = fncForInversed(ans);
return ans;
}
double getDeterm(double x0, double xi, double xik, double p1, double p2)
{
return 1.0 / (xik - xi) * ((x0 - xi) * p2 - (x0 - xik) * p1);
}
double getP(int a, int b, double x0)
{
if (a == b)
return yarr[a];
return getDeterm(x0, xarr[a], xarr[b], getP(a, b - 1, x0), getP(a + 1, b, x0));
}
// многочленом Лагранжа по схеме Эйткена
// Суть метода:
// Критерий сходимости:
double lagrangeApprox(double x0)
{
int n = xarr.Length;
double[,] mtrEik = new double[n, n];
for (int i = 0; i < n; ++i)
{
mtrEik[i, 0] = yarr[i];
mtrEik[i, i] = yarr[i];
}
int nearestA = 0;
int nearestB = 0;
for (int i = 0; i < xarr.Length; ++i)
if (xarr[i] > x0)
{
nearestB = i;
break;
}
nearestA = nearestB - 1;
int pw = 1;
double fDiff = getP(nearestA, nearestA + pw, x0) - getP(nearestA, nearestA, x0);
while (pw + 1 + nearestA < n)
{
double sDiff = getP(nearestA, nearestA + pw + 1, x0) - getP(nearestA, nearestA + pw, x0);
if (Math.Abs(sDiff) > Math.Abs(fDiff))
break;
fDiff = sDiff;
pw++;
}
return getP(nearestA, nearestA + pw, x0);
}
public Form1()
{
InitializeComponent();
}
private void Form1_Load(object sender, EventArgs e)
{
}
private void label2_Click(object sender, EventArgs e)
{
}
private void button1_Click(object sender, EventArgs e)
{
chart1.Series[0].ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Line;
chart1.Series[1].ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Point;
int vl = Convert.ToInt32(textBox4.Text);
// Ньютон
if (vl == 3)
{
getNewtone();
textBox1.Text = ans1.ToString();
textBox2.Text = ans2.ToString();
textBox3.Text = ans3.ToString();
int step = (int)(seriesSize / 1.0 + 1);
step = 1;
for (int i = 0; i < seriesSize; i = i + step)
{
chart1.Series[0].Points.AddXY(dataX[i], dataY[i]);
}
}
// Зейдель
if (vl == 1)
{
getZeidelMethod();
textBox1.Text = ans1.ToString();
textBox2.Text = ans2.ToString();
textBox3.Text = ans3.ToString();
int step = (int)(seriesSize / 1.0 + 1);
step = 1;
for (int i = 0; i < seriesSize; i = i + step)
{
chart1.Series[0].Points.AddXY(dataX[i], dataY[i]);
}
}
// Интерполяция
if (vl == 2)
{
textBox5.Text = findRootByInversedInterpolation(0).ToString();
//findRootByInversedInterpolation(0);
//double step = 0.01;
//for (double a = 0.0; a <= xarr[xarr.Length - 2]; a += step)
//{
// chart1.Series[0].Points.AddXY(a, lagrangeApprox(a));
//}
//chart1.Series[1].MarkerSize = 9;
//for (int i = 0; i < xarr.Length; ++i)
//{
// chart1.Series[1].Points.AddXY(xarr[i], yarr[i]);
//}
}
// Лагранж
if (vl == 4)
{
lagrangeApprox(0);
double step = 0.01;
for (double a = 0.0; a <= xarr[xarr.Length - 2]; a += step)
{
chart1.Series[0].Points.AddXY(a, lagrangeApprox(a));
}
chart1.Series[1].MarkerSize = 9;
for (int i = 0; i < xarr.Length; ++i)
{
chart1.Series[1].Points.AddXY(xarr[i], yarr[i]);
}
}
if (vl == 10)
{
for (int i = 0; i < seriesSize; i++)
{
chart1.Series[0].Points.AddXY(dataX[i], dataY[i]);
}
chart1.Series[1].MarkerSize = 9;
for (int i = 0; i < xarr.Length; ++i)
{
chart1.Series[1].Points.AddXY(xarr[i], yarr[i]);
}
}
}
}
}