Математические модели технических объектов

Также дополнительно  можно установить режим анализа  погрешности, позволяющий проследить изменение погрешности. Результаты по этому анализу заносятся в таблицу на вкладке “Анализ погрешности”:

Для получения результатов необходимо нажать кнопку “Запуск” на вкладке “Решение дифференц. уравнения”. Там же появляются результаты в таблице:

На этой вкладке также  выводится дополнительная информация, такая как

• максимальная погрешность
• количество итераций метода
• число разбиений Nx, Ny и N
• потраченное время

Текст программы

Определение типов данных и методов

 

using System;

using System.IO;

namespace MMKursProject

{

public class Matrix

{

public int N;

public int Nx;

public bool fiveDiag;

public bool econom;

public bool transpon;

public double[] A, B, C;

public double[][] M;

 

public Matrix(int N, int Nx, bool econom, bool fiveDiag)

{

this.Nx = Nx;

this.N=N;

this.fiveDiag=fiveDiag;

this.econom=econom;

transpon=false;

if(econom)

{ A = new double[N];

B = new double[N - 1];

C = new double[N - Nx];

}

else

{ M=new double[N][];

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

M[i]=new double[N];

}

}

public Matrix(int N, int Nx, string str)

{

this.Nx = Nx;

this.N=N;

transpon=false;

if(str=="econom") this.econom=true;

else if(str=="full") this.econom=false;

else throw new ArgumentOutOfRangeException();

if(econom)

{

A = new double[N];

B = new double[N - 1];

C = new double[N - Nx];

}

else

{

M=new double[N][];

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

M[i]=new double[N];

}

}

public Matrix(int N, int Nx, double[] A, double[] B, double[] C)

{

if ((A.Length!=N)||(B.Length!=N-1)||(C.Length!=N-Nx))

throw new ArgumentOutOfRangeException();

this.Nx = Nx;

transpon=false;

this.N=N;

this.econom=true;

this.A = new double[N];

this.B = new double[N - 1];

this.C = new double[N - Nx];

for (int i=0;i<N;i++)  this.A[i]=A[i];

for (int i=0;i<N-1;i++) this.B[i]=B[i];

for (int i=0;i<N-Nx;i++) this.C[i]=C[i];

}

public Matrix(int N, int Nx, Vector A, Vector B, Vector C)

{

if ((A.GetLenght!=N)||(B.GetLenght!=N-1)||(C.GetLenght!=N-Nx))

throw new ArgumentOutOfRangeException();

this.Nx = Nx;

this.N=N;

transpon=false;

this.econom=true;

this.A = new double[N];

this.B = new double[N - 1];

this.C = new double[N - Nx];

for (int i=0;i<N;i++)  this.A[i]=A[i];

for (int i=0;i<N-1;i++) this.B[i]=B[i];

for (int i=0;i<N-Nx;i++) this.C[i]=C[i];

}

public int GetLength

{

get { return A.Length; }

}

public Matrix(StreamReader sr): this(Int16.Parse(sr.ReadLine()),Int16.Parse(sr.ReadLine()),sr.ReadLine())

{ string str;

transpon=false;

int i=0,j=0,k=0,newk=0;

if(!econom)

{ while (((str=sr.ReadLine())!=null)&&(i<N))

{ j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf('\t',k+1);

this[i,j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while(k!=-1);

i++;

}

}

else

{ fiveDiag = true;

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

A[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N));

}

else throw new ArrayTypeMismatchException();

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

B[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N-1));

}

else throw new ArrayTypeMismatchException();

if ((str=sr.ReadLine())!=null)

{

j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

C[j]=double.Parse(newk>0?str.Substring(k,

newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<N-Nx));

}

else throw new ArrayTypeMismatchException();

}

}

 

public double this[int i, int j]

{

get {return GetElement(i, j);}

set {SetElement(i, j, value);}

}

 

public double GetElement(int i, int j)

{ if(transpon)

{ int temp=i;

i=j;

j=temp;

}

if ((i < 0) || (j < 0) || (i >= N) || (j >= N))// Выход за границы

throw new ArgumentOutOfRangeException();

if(econom)

{ if (fiveDiag) // если 5-диагональная матрица

{ if (i==j)// Диагональ

return A[i];

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

{

if(i==j+1) return B[j];

if(i==j+Nx) return C[j];

}

else// Верхний треугольник

{

if(i==j-1) return B[i];

if(i==j-Nx) return C[i];

}

}

else //нижняя треугольная трехдиагональная матрица

{

if (i==j)// Диагональ

return A[i];

// Нижний треугольник

if (i > j)

{

if(i==j+1) return B[j];

if(i==j+Nx) return C[j];

}

 

}

return 0;

}

else return M[i][j];

}

public void SetElement(int i, int j, double value)

{ if(transpon)

{ int temp=i;

i=j;

j=temp;

}

if ((i < 0) || (j < 0) || (i >= N) || (j >= N))// Выход за границы

throw new ArgumentOutOfRangeException();

if(econom)

{ if (fiveDiag)

{ if (i==j)// Диагональ

A[i]=value;

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

{

if(i==j+1) B[j]=value;

if(i==j+Nx) C[j]=value;

}

else// Верхний треугольник

{

if(i==j-1) B[i]=value;

if(i==j-Nx) C[i]=value;

}

}

else

{ if (i==j)// Диагональ

A[i]=value;

// Нижний треугольник

 

if((i==j+1)&&(i>j)) B[j]=value;

if((i==j+Nx)&&(i>j)) C[j]=value;

if((i==j-1)&&(i<j)) B[i]=value;

if((i==j-Nx)&&(i<j)) C[i]=value;

 

}

}

else M[i][j]=value;

}

public static Matrix operator *(Matrix m1, Matrix m2)

// умножение матрицы  на матрицу

{ Matrix res = null;

int len=0;

if ((len=m1.N) != m2.N) throw new ArgumentOutOfRangeException();

else

{ res = new Matrix(len,m1.Nx,false,true);

for (int n = 0; n < len; n++)

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

{ res[n,i]=0;

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

{res[n,i] += m1[n, j]*m2[j,i];}

}

}

}

return res;

}

public override string ToString()

{ string str="";

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

{for(int j=0; j<N; j++) str+=this[i,j].ToString("#,##0.00")+

((j<N-1)?"\t":"");//"#,###0.000"

str+='\n';

}

return str;

}

}

public class Vector

{ public double[] V;

public Vector(int size)

{

V=new double[size];

}

public Vector(int size, double[] V)

{

this.V=new double[size];

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

this.V[i]=V[i];

}

public int GetLenght

{ get { return V.Length; }

}

public Vector(StreamReader sr):this(Int16.Parse(sr.ReadLine()))

{ string str;

int j=0,k=0,newk=0;

if ((str=sr.ReadLine())!=null)

{ j=0;

k=0;

newk=0;

do

{ newk=str.IndexOf(' ',k+1);

V[j]=double.Parse(newk>0?str.Substring(k,newk-k):str.Substring(k));

k=newk;

j++;

}while((k!=-1)&&(j<V.Length));

}

else throw new ArrayTypeMismatchException();

}

 

public double this[int i]

{ get

{ if(i>=V.Length) throw new ArgumentOutOfRangeException();

return V[i];

}

set

{ if(i>=V.Length) throw new ArgumentOutOfRangeException();

V[i]=value;

}

}

public static double operator *(Vector v1, Vector v2)

//умножение вектора  на вектор

{ double res = 0;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

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

{ res += v1[i]*v2[i];

}

return res;

}

public static Vector operator +(Vector v1, Vector v2)

// сумма двух векторов

{ Vector res = null;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

{ res = new Vector(len);

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

{res[i] = v1[i] + v2[i];

}

}

return res;

}

public static Vector operator -(Vector v1, Vector v2)

// разность двух  векторов

{

Vector res = null;

int len=0;

if ((len=v1.V.Length) != v2.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

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

{

res[i] = v1[i] - v2[i];

}

}

return res;

}

public static Vector operator *(Matrix m, Vector v)

// умножение матрицы  на вектор

{

Vector res = null;

int len=0;

if ((len=m.N) != v.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

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

{

res[i] = 0;

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

{res[i] += m[i, j]*v[j];

}

}

}

return res;

}

public static Vector operator *(Vector v, Matrix m)

// умножение вектора  на матрицу

{ Vector res = null;

int len=0;

if ((len=m.N) != v.V.Length) throw new ArgumentOutOfRangeException();

else

{

res = new Vector(len);

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

{ res[j] = 0;

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

{res[j] += m[i, j]*v[i];

}

}

}

return res;

}

public static Vector operator *(double d, Vector v)

//умножение скаляра  на вектор

{ int len=v.V.Length;

Vector res = new Vector(len);

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

{res[i] = v[i]*d;

}

return res;

}

public static Vector operator *(Vector v, double d)

//умножение вектора  на скаляр

{ int len=v.V.Length;

Vector res = new Vector(len);

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

{res[i] = v[i]*d;

}

return res;

}

public static double Abs(Vector v)

{return (Math.Sqrt((v*v)));

}

 

public override string ToString()

{ string str="";

int len=V.Length;

for(int i=0;i<len;i++) str+=V[i].ToString("#,#####0.00000")+" ";

return str;

}

}

}

 

Метод сопряженных градиентов

 

using System;

using System.Windows.Forms;

namespace MMKursProject

{

public class Solving

{ public static Vector X;

public static int Nx = 4;

public static int Ny = 4;

public static double Eps=0.00000001;

public static double xLeft=0, xRight=1, yLeft=0, yRight=1;

public static double hx, hy;

public static int N,i,j;

public static double kappa1 = 1;

public static double kappa2 = 1;

public static double kappa3 = 1;

public static double kappa4 = 1;

public static double[] xi, yj;

public static double[] Ck, Bk, Ak, Fk,Uk,XX;

public static int k = 0, iter = 0;

public static double maxerror = 0;

public static DateTime start;

public static TimeSpan ts;

public static int State=0;

 

public static void Solve()

{

try

{

start = DateTime.Now;

N = Nx*Ny;

hx = (xRight - xLeft)/Nx;

hy = (yRight - yLeft)/Ny;

xi = new double[Nx + 1];

yj = new double[Ny + 1];

// сетка по х

for (i = 0; i <= Nx; i++)

{ xi[i] = xLeft + i*hx;  }

// сетка по у

for (j = 0; j <= Ny; j++)

{ yj[j] = yLeft + j*hy;  }

Ck = new double[N];

Bk = new double[N];

Ak = new double[N];

Fk = new double[N];

Uk = new double[N];

XX = new double[N];

for (j=0; j<=Ny-1; j++)

for (i=0; i<=Nx-1; i++)

{

k=j*Nx+i;

Uk[k]=u(xi[i],yj[j]);

}

 

Ak[0]=(hy/2)*(k1(xi[0],xi[1],yj[0])/hx + kappa1) + (hx/2)*(k2(xi[0],yj[0],yj[1])/hy + kappa3);

Fk[0]=(hx*hy/4)*fun(xi[0],yj[0]) + (hy/2)*V1(yj[0]) + (hx/2)*V3(xi[0]);

for ( i = 1; i<=(Nx-1); i++)

for (j = 1; j<= (Ny-1); j++)

{

k = j*Nx + i;

Ak[k]= (hy/hx)*(k1(xi[i],xi[i+1],yj[j]) + k1(xi[i],xi[i-1],yj[j])) + (hx/hy)*(k2(xi[i],yj[j],yj[j+1]) + k2(xi[i],yj[j],yj[j-1]));

Bk[k]= -(hy/hx)*k1(xi[i],xi[i-1],yj[j]) ;

Ck[k]= -(hx/hy)*k2(xi[i],yj[j],yj[j-1]);

Fk[k]= hx*hy*fun(xi[i],yj[j]);

if (i==(Nx-1)) Fk[k]= Fk[k] + (hy/hx)*V2(yj[j])*k1(xi[i],xi[i+1],yj[j]);

if (j==(Ny-1)) Fk[k]= Fk[k] + (hx/hy)*V4(xi[i])*k2(xi[i],yj[j],yj[j+1]); 

}

for (j = 1; j<=(Ny-1); j++)

{

i = 0;

k=j*Nx + i;

Ak[k]=hy*(k1(xi[i],xi[i+1],yj[j])/hx + kappa1) + (hx/(2*hy))*(k2(xi[i],yj[j],yj[j+1]) + k2(xi[i],yj[j],yj[j-1]));

Ck[k]= -(hx/(2*hy))*k2(xi[i],yj[j],yj[j-1]);

Fk[k]=(hx/2)*hy*fun(xi[i],yj[j]) + hy*V1(yj[j]);

if (j==(Ny-1)) Fk[k]= Fk[k] + (hx/(2*hy))*V4(xi[i])*k2(xi[i],yj[j],yj[j+1]);

}

for (i = 1; i<=(Nx-1); i++)

{

j = 0;

k=j*Nx + i;

Ak[k]=hx*(k2(xi[i],yj[j],yj[j+1])/hy + kappa3) + (hy/(2*hx))*(k1(xi[i],xi[i+1],yj[j]) + k1(xi[i],xi[i-1],yj[j]));

Bk[k]= -(hy/(2*hx))*k1(xi[i],xi[i-1],yj[j]);

Fk[k]=hx*(hy/2)*fun(xi[i],yj[j]) + hx*V3(xi[i]);

if (i==(Nx-1)) Fk[k]= Fk[k] + (hy/(hx*2))*V2(yj[j])*k1(xi[i],xi[i+1],yj[j]);

}

 

//------------------------------------------------------------------------------------

// формируем вектора  для матрицы

Vector A = new Vector(N);

Vector B = new Vector(N - 1);

Vector C = new Vector(N - Nx);

Vector F = new Vector(N);

Vector X0 = new Vector(N);

Vector UU = new Vector(N);

for (k = 0; k < N; k++) A[k] = Ak[k];

for (k = 0; k < N - 1; k++) B[k] = Bk[k + 1];

for (k = 0; k < N - Nx; k++) C[k] = Ck[k + Nx];

for (k = 0; k < N; k++) F[k] = Fk[k];

for (k = 0; k < N; k++) UU[k] = Uk[k]; 

// формируем матрицу

Matrix MM = new Matrix(N, Nx, A, B, C);

MM.fiveDiag = true;

//X = Metod(MM, F, X0); // метод сопряженных градиентов без предобусловливания

X = MetodHoleckogo(MM, F, X0); // метод сопряженных градиентов без предобусловливания

maxerror=Math.Abs(UU[0]-X[0]);

for (i=1;i<N;i++)

{

if(maxerror<Math.Abs(UU[i]-X[i]))

{ maxerror=Math.Abs(UU[i]-X[i]); }

}

DateTime end = DateTime.Now;

ts=end-start;

State=2;

 

}

catch (Exception exc)

{ MessageBox.Show(exc.Message+"***");  }

}

// метод сопряженных градиентов без предобуславливания

public static Vector Metod(Matrix A, Vector F, Vector x0)

{ k = 0;

bool Stop = false;

Vector r0 = new Vector(F.GetLenght);

Vector r1 = new Vector(F.GetLenght);

Vector S1 = new Vector(F.GetLenght);

Vector x = new Vector(F.GetLenght);

Vector x00 = new Vector(F.GetLenght);

Double alfa, betta;

x00 = x0;

k = 0;

// невязка

r0 = (F - A*x0);

//  вспомогательный  вектор в методе

S1 = r0;

while (!Stop)

{

k++;

// подсчет альфы

alfa = (r0*r0)/((A*S1)*S1);

x = x0 + alfa*S1;

r1 = r0 - alfa*(A*S1);

// проверка на  сходимость - завершение метода

if ((Vector.Abs(r1)/Vector.Abs(F)) < 0.000001)

{ if ((Vector.Abs(F - A*x00)/Vector.Abs(F)) > 0.000001)

{ Stop = true;  }

else return x;

}

betta = (r1*r1)/(r0*r0);

S1 = r1 + betta*S1;

r0 = r1;

x0 = x;

}

return x;

}

// метод сопряженных градиентов с предобуславливанием

public static Vector MetodHoleckogo(Matrix A, Vector F, Vector x0)

{

iter = 0;

bool Stop = false;

Vector r0 = new Vector(F.GetLenght);

Vector r1 = new Vector(F.GetLenght);

Vector S1 = new Vector(F.GetLenght);

Vector x = new Vector(F.GetLenght);

Vector x00 = new Vector(F.GetLenght);

Vector w0 = new Vector(F.GetLenght);

Vector w1 = new Vector(F.GetLenght);

Matrix L;

double alfa, betta;

x00 = x0;

iter = 0;

// невязка

r0 = (F - A*x0);

L = CreateMatrixL(A);

w0 = FindW(L, r0);

//  вспомогательный  вектор в методе

S1 = w0;

while (!Stop)

{

iter++;

// подсчет альфы

alfa = (w0*r0)/((A*S1)*S1);

x = x0 + alfa*S1;

r1 = r0 - alfa*(A*S1);

// проверка на  сходимость - завершение метода

if ((Vector.Abs(r1)/Vector.Abs(F)) < 0.000001)

{

if ((Vector.Abs(F - A*x00)/Vector.Abs(F)) > 0.000001)

{Stop = true;}

else return x;

}

w1 = FindW(L, r1);

betta = (w1*r1)/(w0*r0);

S1 = w1 + betta*S1;

r0 = r1;

w0 = w1;

x0 = x;

}

return x;

}

public static Matrix CreateMatrixL(Matrix M)

{ //Получаем вектора a,b,c из основной матрицы для формирования векторов в матрицу L

Vector A = new Vector(M.GetLength);

Vector B = new Vector(M.GetLength);

Vector C = new Vector(M.GetLength);

double[] Al = new double[M.GetLength];

double[] Bl = new double[M.GetLength];

double[] Cl = new double[M.GetLength];

for (int i = 0; i < M.N; i++)

{ for (int j = 0; j < M.N; j++)

{ if (i == j) // Диагональ

A[i] = M[i, j];

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

{

if (i == j + 1) B[j] = M[i, j];

if (i == j + M.Nx) C[j] = M[i, j];

}

}

}

Al[0]=Math.Sqrt(A[0]);

Bl[0]=B[0]/Al[0];

Cl[0] = C[0]/Al[0];

for (int i = 1; i < M.Nx; i++)

{ Al[i] = Math.Sqrt(A[i]-Bl[i-1]*Bl[i-1]);

Bl[i] = B[i]/Al[i];

Cl[i] = C[i]/Al[i];

}

for (int i = M.Nx; i < M.N; i++)

{ Al[i] = Math.Sqrt(A[i]-Bl[i-1]*Bl[i-1]-Cl[i-M.Nx]*Cl[i-M.Nx]);

Bl[i] = B[i]/Al[i];

Cl[i] = C[i]/Al[i];

}

//формируем вектора  для построения матрицы

A = new Vector(M.N);

B = new Vector(M.N - 1);

C = new Vector(M.N - M.Nx);

for (int i = 0; i < M.N; i++) A[i] = Al[i];

for (int i = 0; i < M.N - 1; i++) B[i] = Bl[i];

for (int i = 0; i < M.N - M.Nx; i++) C[i] = Cl[i];

Matrix L = new Matrix(M.N, M.Nx, A, B, C);

return L;

}

// метод  Гаусса - решает систему вниз

public static Vector GaussDown(Matrix M, Vector F)

{ Vector x = new Vector(M.N);

x[0] = F[0]/M.GetElement(0, 0);

for (int i = 1; i < M.Nx; i++) x[i] = (F[i] - M.GetElement(i, i - 1)*x[i - 1])/M.GetElement(i, i);

for (int i = M.Nx; i < M.N; i++)

x[i] = (F[i]-M.GetElement(i, i-M.Nx)*x[i-M.Nx]-M.GetElement(i,i-1)*x[i-1])/M.GetElement(i, i);

return x;

}

// метод  Гаусса - решает систему вверх

public static Vector GaussUp(Matrix M, Vector F)

{ M.transpon = true;

Vector x = new Vector(M.N);

x[M.N - 1] = F[M.N - 1]/M.GetElement(M.N - 1, M.N - 1);

for (int i = M.N - 2; i >= (M.N - M.Nx); i--)

x[i] = (F[i] - M.GetElement(i, i + 1)*x[i + 1])/M.GetElement(i, i);

for (int i = M.N - M.Nx - 1; i >= 0; i--)

x[i] = (F[i]-M.GetElement(i,i+M.Nx)*x[i+M.Nx]-M.GetElement(i,i+1)*x[i+1])/M.GetElement(i,i);

M.transpon = false;

return x;

}

// нахождение дополнительного вектора w

public static Vector FindW(Matrix L, Vector r)

{ Vector temp = new Vector(L.GetLength, GaussDown(L, r).V);

return GaussUp(L, temp);

}

// функция к1

public static double k1(double xi, double xi_1, double yi)

{ double x = (xi + xi_1)/2;

double res;

//res = 2; //test 1

//res = 2*x + yi + 1; //test 2

res = 1 + 2*x*x + yi; //test3

return res;

}

// функция к2

public static double k2(double xi, double yi, double yi_1)

{ double res;

double y = (yi + yi_1)/2;

//res = 5; //test 1

//res = xi + 2*y + 1; //test2;

res = 1 + xi + 2*y*y; //test3;

return res;

}

//решение уравнения - u

public static double u(double x, double y)

{

double res;

//res = 3; // test1

//res = 2*x+3*y+5; //test2

res = 2*x*x + 3*y*y*y + 4; //test3;

return res;

}

 

public static double V1(double y)

{

double res;

//res = 3; // test1

//res =kappa1*(2*xLeft+3*y+5) - (2*xLeft+y+1)*2; //test2

res = kappa1*(2*xLeft*xLeft+3*y*y*y+4)-(1+2*xLeft*xLeft+y)*(4*xLeft);//test3

return res;

}

 

public static double V2(double y)

{

double res;