Основы компьютерного моделирования физических и технических систем. Математическую модель указанной физической системы, страница 3

string footor = "</Table> \n <WorksheetOptions xmlns=\"urn:schemas-microsoft-com:office:excel\"> \n  <PageSetup>  \n  <Header x:Margin=\"0.3\"/>   \n <Footer x:Margin=\"0.3\"/>    \n<PageMargins x:Bottom=\"0.75\" x:Left=\"0.7\" x:Right=\"0.7\" x:Top=\"0.75\"/>   \n</PageSetup>   \n<Unsynced/>   \n<Selected/>   \n<Panes>    \n<Pane>     \n<Number>3</Number>     \n<ActiveRow>2</ActiveRow>     \n<ActiveCol>3</ActiveCol>    \n</Pane>   \n</Panes>   \n<ProtectObjects>False</ProtectObjects>   \n<ProtectScenarios>False</ProtectScenarios>  \n</WorksheetOptions> \n</Worksheet> \n</Workbook> ";

StreamWriter Fil;

//удалить

Fil = new StreamWriter("T.xml");

Fil.WriteLine(heder);

Fil.WriteLine("<Table ss:ExpandedColumnCount=\"{0}\" ss:ExpandedRowCount=\"{1}\" x:FullColumns=\"1\"   \nx:FullRows=\"1\" ss:DefaultRowHeight=\"15\">", u.GetLength(0) + 1, u.GetLength(1) + 1);

for (int i = 0; i < u.GetLength(0); i++)

{

Fil.WriteLine("<Row ss:AutoFitHeight=\"0\">\n");

for (int j = 0; j < u.GetLength(1); j++)

{

Fil.Write("<Cell><Data ss:Type=\"Number\">");

string ss = "";

ss = string.Format("{0,9:f14}", u[i, j]);

ss = ss.Replace(',', '.');

Fil.Write(ss);

Fil.WriteLine("</Data></Cell>");

}

Fil.WriteLine("</Row>\n");

}

Fil.WriteLine(footor);

Fil.Close();

long start = DateTime.Now.Ticks;

// Lm = new double[u.GetLength(0), u.GetLength(1)];

// D = new double[u.GetLength(0), u.GetLength(1)];

StreamWriter ffii = new StreamWriter("A.txt");

StreamWriter ffi2 = new StreamWriter("B.txt");

for (int i = 0; i < u.GetLength(0); i++)

{

for (int j = 0; j < u.GetLength(1); j++)

{

ffii.Write(u[i, j] + " ");

}

ffii.WriteLine();

ffi2.WriteLine(b[i]);

}

ffii.Close();

ffi2.Close();

//Решение СЛАУ

mL = new double[u.GetLength(0), u.GetLength(1)];

mD = new double[u.GetLength(0), u.GetLength(1)];

mU = new double[u.GetLength(0), u.GetLength(1)];

for (int i = 0; i < u.GetLength(0); i++)

{

mL[i,i] = 1;

mU[i,i] = 1;

}

//находим первый столбец L[][] и первую строку D[][]

mD[0,0] = u[0,0];

for (int i = 0; i < u.GetLength(0); i++)

{

mL[i,0] = u[i,0] / mD[0,0];

mU[0,i] = u[0,i] / (mL[0,0] * mD[0,0]);

}

//дальше вычисляем L[][], D[][] по формуле

double sum = 0;

for (int i = 1; i < u.GetLength(0); i++)

{

for (int j = 1; j < u.GetLength(0); j++)

{

if (i == j) //нижний треугольник

{

sum = 0;

for (int k = 0; k < j; k++)

sum += mL[i,k] * mU[k,j] * mD[k,k];

mD[i,i] = (u[i,i] - sum);///D[i][i];

}

else

if (i > j)

{

sum = 0;

for (int k = 0; k < i; k++)

sum += mL[i,k] * mU[k,j] * mD[k,k];

mL[i,j] = (u[i,j] - sum) / mD[j,j];

}

else// верхний

{

sum = 0;

for (int k = 0; k < i; k++)

sum += mL[i,k] * mU[k,j] * mD[k,k];

mU[i,j] = (u[i,j] - sum) / (mD[i,i] * mL[i,i]);

}

}

}

int ln = u.GetLength(0);  //размерность матрицы

double[] vY = new double[ln];

double[] vZ = new double[ln];

vX = new double[ln]; //инициализация массива ответов

double hh=0;

for (int i = 0; i <ln; i++)

{

hh = b[i];// hh присваиваем ответ из массива сил

for (int j = 0; j <= i; j++)

hh = hh - vY[j] * mL[i, j]; //вычетание элементов строки матрицы

vY[i] = hh / mL[i, i];//hh  делим на коэффициент при иксу

}

for (int i = 0; i < ln; i++)

{

//hh = b[i];// hh присваиваем ответ из массива сил

//for (int j = 0; j < ln; j++)

//hh = hh - vZ[j] * mD[i, j]; //вычетание элементов строки матрицы

vZ[i] = vY[i] / mD[i, i];//hh  делим на коэффициент при иксу