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 NumMethods3
{
public partial class Form1 : Form
{
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];
}
}
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;
}
public Form1()
{
InitializeComponent();
}
double[] xval;
double[] yval;
void limDifferencesMethod()
{
double a = 1.2;
double b = 1.9;
double h = 0.05;
double A = 5.1;
double B = 1.4;
double al1 = 2;
double al2 = -1.7;
double bet1 = -4;
double bet2 = 5;
double p = -2;
double q = -3;
int n = (int) ((b - a) / h);
xval = new double[n + 1];
yval = new double[n + 1];
Mtr matrix = new Mtr(n + 1, n + 1);
Mtr fun = new Mtr(n + 1, 1);
for (int i = 1; i < n; ++i)
{
matrix[i, i - 1] = (1 - h / 2.0 * p);
matrix[i, i] = -(2 - h * h * q);
matrix[i, i + 1] = (1 + h / 2.0 * p);
fun[i, 0] = h * h * Math.Sin(3 * (a + h * i));
}
fun[0, 0] = 2 * A * h;
fun[n, 0] = 2 * B * h;
matrix[0, 0] = 2 * h * al1 - 3 * al2;
matrix[0, 1] = 4 * al2;
matrix[0, 2] = -al2;
matrix[n, n - 2] = bet2;
matrix[n, n - 1] = - 4 * bet2;
matrix[n, n] = (2 * h * bet1 + 3 * bet2);
Mtr ans = inversedLuSolver(matrix, fun);
for (int i = 0; i <= n; ++i)
{
xval[i] = a + h * i;
yval[i] = ans[i, 0];
}
}
double funcForRunge(double x, double y)
{
return x * y;
}
double rngP(double x)
{
return -2;
}
double rngQ(double x)
{
return -3;
}
double rngFx(double x)
{
return Math.Sin(3 * x);
}
double[] valuesDelta = new double[1000];
double[] valuesGamma = new double[1000];
double[] valuesPassAnswer = new double[1000];
double funcDelta(double x, double y, int ind)
{
return -rngP(x) * y - y * y - rngQ(x);
}
double funcGamma(double x, double y, int ind)
{
return -(rngP(x) + valuesDelta[ind]) * y + rngFx(x);
}
double finalFncForRunge(double x, double y, int ind)
{
return valuesDelta[ind] * y + valuesGamma[ind];
}
delegate double fnc(double x, double y, int ind);
double runge_kutta(double b, double y0, double x0, double h, fnc funcForRunge, double[] output, bool reverse = false)
{
double x = x0;
double y = y0;
int ind = (int)((b - x0) / h);
int dind = -1;
if (!reverse)
{
ind = 0;
dind = 1;
}
if (!reverse)
{
while (x < b)
{
double ni1 = funcForRunge(x, y, ind);
double ni2 = funcForRunge(x + h / 2.0, y + h / 2.0 * ni1, ind);
double ni3 = funcForRunge(x + h / 2.0, y + h / 2.0 * ni2, ind);
double ni4 = funcForRunge(x + h, y + h * ni3, ind);
double delta = h / 6.0 * (ni1 + 2 * ni2 + 2 * ni3 + ni4);
y = y + delta;
x = x + h;
output[ind] = y;
ind += dind;
}
}
else
{
while (b < x)
{
double ni1 = funcForRunge(x, y, ind);
double ni2 = funcForRunge(x + h / 2.0, y + h / 2.0 * ni1, ind);
double ni3 = funcForRunge(x + h / 2.0, y + h / 2.0 * ni2, ind);
double ni4 = funcForRunge(x + h, y + h * ni3, ind);
double delta = h / 6.0 * (ni1 + 2 * ni2 + 2 * ni3 + ni4);
y = y + delta;
x = x + h;
output[ind] = y;
ind += dind;
}
}
return y;
}
double methodDiffPass()
{
double a = 1.2;
double b = 1.9;
double h = 0.05;
int steps = (int)((b - a) / h);
double A = 5.1;
double B = 1.4;
double al1 = 2;
double al2 = -1.7;
double bet1 = -4;
double bet2 = 5;
runge_kutta(b, -al1 / al2, a, h, funcDelta, valuesDelta);
runge_kutta(b, A / al2, a, h, funcGamma, valuesGamma);
runge_kutta(a, (B - bet2 * valuesGamma[steps - 1]) / (bet1 + bet2 * valuesDelta[steps - 1]) , b, -h, finalFncForRunge, valuesPassAnswer, true);
return 0;
}
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 getByHoin(double a, double b, double h, fnc function, double y0)
{
double y = y0;
while (a < b)
{
double fnc = function(a, y, 0);
y = y + h / 2.0 * (fnc + function(a + h, y + h * fnc, 0));
a += h;
}
return y;
}
double[] adamsOrd = new double[1000];
double[] adamsAbc = new double[1000];
double[] adamsFunc = new double[1000];
double shootDeriv(double x, double y, int ind)
{
return adamsOrd[ind];
}
double funcForStrelb(double x, double y, int ind)
{
return Math.Exp(-x) - y;
}
/* double adamsMethod(double a, double b, double h, double y0, double derv)
{
double y = y0;
double cpos = a + h;
adamsAbc[0] = a;
adamsOrd[0] = derv;
adamsFunc[0] = funcForStrelb(a, y0, 0);
int i = 0;
for (i = 1; i <= 4; ++i)
{
adamsAbc[i] = cpos;
adamsOrd[i] = getByHoin(a, cpos, h, functForSterlb, y0);
adamsFunc[i] = functForSterlb(adamsAbc[i], adamsOrd[i], 0);
cpos += h;
}
i--;
while (cpos < b)
{
y = 2 * adamsOrd[i] - adamsOrd[i - 1] + h * h * (adamsFunc[i] + 1.0 / 12.0 * getDelta(2, i - 2, adamsFunc) + 1.0 / 12.0 * getDelta(3, i - 3, adamsFunc)
+ 19.0 / 240.0 * getDelta(4, i - 4, adamsFunc));
i++;
cpos += h;
adamsAbc[i] = cpos;
adamsOrd[i] = y;
adamsFunc[i] = functForSterlb(adamsAbc[i], adamsOrd[i], 0);
}
return adamsOrd[i];
}
void shootMethod()
{
}
*/
private void Form1_Load(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;
chart1.Series[2].ChartType = System.Windows.Forms.DataVisualization.Charting.SeriesChartType.Line;
limDifferencesMethod();
methodDiffPass();
for (int i = 0; i < xval.Length; ++i)
{
chart1.Series[0].Points.AddXY(xval[i], yval[i]);
double x = xval[i];
double vl = 5.20377 * Math.Exp(-x) + 0.00280526 * Math.Exp(3 * x) +0.0333333 * Math.Cos(3 * x) - 0.0666667 * Math.Sin(3 * x);
// double vl = Math.Exp(-x) * (4.82068 + 0.00239973 * Math.Exp(4 * x));
chart1.Series[1].Points.AddXY(xval[i], vl);
chart1.Series[2].Points.AddXY(x, valuesPassAnswer[i]);
}
}
}
}