Решение эллиптических уравнений

Министерство образования и науки РФ

Новосибирский Государственный Технический  Университет

 

 

 

 

 

Кафедра Прикладной Математики

 

 

 

 

 

 

Курсовой проект

по дисциплине «Численные методы» 

 

 

 

 

Факультет: ПМИ

Группа:

Студентка: 

Преподаватель:

Вариант : 25

 

 

 

   

 

 

 

 

 

 

 

 

 

 

 

Новосибирск 2010

 

Задание

МКЭ для двумерной краевой  задачи для эллиптического уравнения  в декартовой системе координат. Базисные функции линейные на треугольниках. Краевые условия всех типов. Коэффициент  разложить по линейным базисным функциям. Матрицу СЛАУ генерировать в разреженном строчном формате. Для решения СЛАУ использовать МСГ или ЛОС с неполной факторизацией.

Вариационная постановка

Эллиптическая краевая задача для функции u определяется дифференциальным уравнением

 

заданным в некоторой  области  с границей и краевыми условиями

 

 

 

в которых – значение искомой функции u на границе S_i , а – значение на S_i производной функции u по направлению внешней нормали к поверхности S_i .

Дифференциальное уравнение   (**)  для двумерной эллиптической  краевой задачи в декартовой системе  координат {х.у} может быть записана в виде : 
- + = . 

Базисные функции линейные на треугольниках.

На каждом конечном элементе определим три локальные базисные функции

.

такие. что  равна единице в вершине с координатами (х1 .у1) и нулю в двух других вершинах. Определенные таким образом локальные базисные функции на треугольнике фактически являются - координатами этого треугольника. Поэтому коэффициенты могут быть вычислены по формуле:

 

 

.

Для вычисления интегралов будем пользоваться формулой:

 

Локальная матрица жесткости

  (αⁱ₁α^j_{1 +} αⁱ₂α^j₂ )* |det D|/2

Локальная матрица массы

 

Разложим  по базисным функциям:

=

 .

Вектор правой части

 

 

 

Сборка  глобальной матрицы и глобального  вектора правой части.

При формировании глобальной матрицы из локальных. полученных суммированием соответствующих матриц массы и жесткости, учитываем соответствие локальной и глобальной нумераций каждого конечного элемента. Глобальная нумерация каждого конечного элемента однозначно определяет позиции вклада его локальной матрицы в глобальную. Поэтому, зная глобальные номера соответствующих узлов конечного элемента, определяем и то, какие элементы глобальной матрицы изменятся при учете текущего конечного элемента. Аналогичным образом определяется вклад локального вектора правой части в глобальный. При учете текущего локального вектора изменятся те элементы глобального вектора правой части, номера которых совпадают с глобальными номерами узлов. присутствующих в этом конечном элементе.

Учет первых краевых условий

 

Вместо того. чтобы на место фиктивного i-го уравнения ставить уравнение , диагональный элемент a_ii в матрице заменим некоторым большим числом С, а i-ую компоненту вектора правой части заменим значение С.

Учет вторых краевых условий

 

Данное краевое условие  вносит вклад только в правую часть  СЛАУ, т.е.в соответствии с тем, на каком ребре задано второе краевое условие, в правую часть к соответствующим номерам добавляется

B^S2 =   * , где – длина ребра ((x1.y1) и (x2.y2) – координаты вершин ребра )

Учет третьих краевых условий

 

Третье краевое условие  дает вклад и в матрицу и  в правую часть СЛАУ:

 

Тогда  , где длина ребра. на котором задано краевое условие

 

Решение СЛАУ

Для решения СЛАУ с полученной матрицей в разреженно-строчном формате  будем применять локально-оптимальную схему с предобусловливанием неполным разложением Холесского.

Представление исходных данных

* В файле nvtr.txt . хранится  количество конечных элементов. областей. узлов сетки и глобальная нумерация узлов в порядке локальной

           * В файле nvkat.txt . хранится номер подобласти которой принадлежит i-тый конечный элемент.

           * В файле koord.txt . хранятся координаты узлов в порядке глобальной нумерации

*В файлах ig.txt и jg.txt хранится  портрет глобальной матрицы в  разреженном строчном формате.

 

 

 

 

 

 

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

 

#include <stdio.h>

#include <conio.h>

#include <math.h>

#define N 10

#define k_max 100000 // максимальное количество итераций

#define epsilon 1e-40 // точность

const double C=1e+20; // большое-большое число

double xy[N][2], alfa_mat[3][3], local_gest[3][3], local_mas[3][3], local[3][3], local_f[3], *glob_f, *gg, *diag;

int kol_el, kol_obl, kol_uzl,  nvtr[N][3], nvkat[N], *ig, *jg;//кол-во кон.эл-тов, областей, узлов

double *x,*r,*p,*q,*z,*diag_,*gg_,*f;

// Подпрограмма зануления  вектора х

void x_null()

{

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

x[i]=0.0;

}

// Подпрограмма перемножения  матрицы на вектор

void multi(double *f,double *x)

{

int i,j,l,i0;

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

{

f[i]=diag[i]*x[i];

}

 

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

{

int kolvo=ig[i+1]-ig[i];

for(l=0,i0=ig[i] ; l<kolvo ; l++,i0++)

{

f[i]+=gg[i0]*x[jg[i0]];

f[jg[i0]]+=gg[i0]*x[i];

}

}

}

// Подпрограмма подсчета  скалярного произведения

double skalyar(double *x, double *y)

{

double rezult=0.0;

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

{

rezult+=x[i]*y[i];

}

return rezult;

}

// Подпрограмма неполного  разложения Холесского (LL_T)

void l_mat()

{

int i,j,jk,ik,l;

double sum,sum_diag;

 

for (i=0;i<kol_uzl;i++) // Переприсваивание элементов матрицы gg

{

diag_[i]=diag[i];

}

for (i=0;i<ig[kol_uzl];i++)

{

gg_[i] = gg[i];

}

 

for (i=0;i<kol_uzl;i++) // Само LL_T

{

sum_diag = 0;

for (j=ig[i] ; j<ig[i+1] ; j++) //Lij , Lji

{

sum = 0; 

jk = ig[jg[j]]; // Над текущим элементом (над диагональю)

ik = ig[i]; // До текущего (на той же строке, слева)

while ( (ik<j) && ( jk< ig[jg[j]+1] ) ) // Пока на дошли дотекущего в этой строке или конца верхнего профиля

{

l=jg[jk]-jg[ik]; //Сравнить положения элементов в строке и "столбце"

if (l==0) //Если совпали, то перемножаем и плюсуем

{

sum+= gg_[jk] * gg_[ik];

ik++;

jk++;

}

jk += (l<0); //Иначе передвинуться в сторону большего индекса

ik += (l>0);

}

gg_[j]-=sum;

gg_[j]/=diag_[jg[j]]; // Окончательный подсчет для внедиаг-го эл-та

sum_diag+= gg_[j] * gg_[j];

}

diag_[i]-=sum_diag;

diag_[i]=sqrt(diag_[i]); // Окончательный подсчет для диаг-го эл-та

}

}

 

// Подпрограмма прямого хода

void pryamoj_hod(double *r)

{

int i,j;

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

{

for (j=ig[i]; j<ig[i+1]; j++)

{

r[i]-= gg_[j] * r[jg[j]];

}

r[i]=r[i]/diag_[i];

}

}

 

// Подпрограмма обратного хода

void obratnyj_hod(double *r)

{

int i,j;

for (j=kol_uzl-1;j>=0;j--)

{

r[j]=r[j]/diag_[j];

for (i=ig[j]; i<ig[j+1]; i++)

{

r[jg[i]]-= gg_[i] * r[j];

}

}

}

// Подпрограмма ЛОС с  предобусловливанием неполным разложением  Холесского

void los_hol(double *x, double *f)

{

//unsigned int time=GetTickCount();

FILE *fp3;

fp3 = fopen("Los_hol_Test.txt","w");

double alf,bet,kvadr_nevyazka=epsilon+1.;

int i,k;

 

x_null(); // Начальные значения

pryamoj_hod(f);

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

{

r[i]=f[i];

}

obratnyj_hod(f);

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

{

z[i]=f[i];

}

multi(p,z);

pryamoj_hod(p);

 

for(k=1; k < k_max && epsilon < kvadr_nevyazka ;k++) // Все самое интересное

{

double skp=skalyar(p,p);

alf=skalyar(p,r) / skp;  

kvadr_nevyazka=skalyar(r,r) - alf*alf*skp;

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

{

x[i]+=alf*z[i];

r[i]-=alf*p[i];

}

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

{

f[i]=r[i];

}

obratnyj_hod(f);

multi(q,f);

pryamoj_hod(q);

bet = -skalyar(p,q) / skp;

 

printf("\n");

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

{

z[i]=f[i] + bet*z[i];

p[i]=q[i] + bet*p[i];

}

 

}

 

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

{

 

fprintf(fp3,"%.15e\n",x[i]);

}

fclose(fp3);

}

//ввод данных

void vvod()

{

FILE *f,*f1;

f=fopen("nvtr.txt","r");//хранит кол-во кон.эл-тов, областей, узлов сетки и нумерацию глобальных узлов в порядке локальной

fscanf(f,"%d %d %d", &kol_el, &kol_obl, &kol_uzl);

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

{

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

{

fscanf(f,"%d",&nvtr[i][j]);

}

}

 

f1=fopen("nvkat.txt","r");//какой подобласти пренадлежит i-тый кон.элемент

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

{

fscanf(f1,"%d",&nvkat[i]);

}

fclose(f);

fclose(f1); 

}

 

//считывание координат  х, у

void vvod_koord()

{

FILE *f2;

f2=fopen("koord.txt","r");

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

{

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

{

fscanf(f2,"%lf",&xy[i][j]);

}

}

fclose(f2);

}

// Аналитика

void analitika(double *u)

{

FILE *f;

f=fopen("u_an.txt","w");

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

{

u[i]=xy[i][0]*xy[i][0]; 

fprintf(f,"%.15e\n",u[i]);  //аналитическая функция

}

fclose(f);

}

//лямбда

int lyam()

{

return 1;

}

 

//гамма

int gamma(int number_obl)

{

return 1;

 

}

//коэффициенты ф

double koef_f(double x,double y)

{

return x*x-2;   

}

//зануление

void null_all()

{

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

{

local_f[i]=0.0; 

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

{

local_gest[i][j]=local_mas[i][j]=local[i][j]=0.0;

}

}

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

{

diag[i]=0.0;

glob_f[i]=0.0; 

}

 

for(int i=0;i<ig[kol_uzl];i++)

gg[i]=0.0;

 

}

//Определитель

double det_D(int i)

{

double detD;

detD=(xy[nvtr[i][1]][0]-xy[nvtr[i][0]][0])*(xy[nvtr[i][2]][1]-xy[nvtr[i][0]][1])-(xy[nvtr[i][2]][0]-xy[nvtr[i][0]][0])*(xy[nvtr[i][1]][1]-xy[nvtr[i][0]][1]);

return detD;

}

//длина ребра

double dlina(int uzl1, int uzl2)

{

double l=(xy[uzl2][0]-xy[uzl1][0])*(xy[uzl2][0]-xy[uzl1][0])+(xy[uzl2][1]-xy[uzl1][1])*(xy[uzl2][1]-xy[uzl1][1]);

l=sqrt(l);

return l;

}

//Коэффициенты альфа

void alfa(int i)

{

double w=det_D(i);

alfa_mat[0][0]=(xy[nvtr[i][1]][0]*xy[nvtr[i][2]][1]-xy[nvtr[i][2]][0]*xy[nvtr[i][1]][1])/w;

alfa_mat[0][1]=(xy[nvtr[i][1]][1]-xy[nvtr[i][2]][1])/w;

alfa_mat[0][2]=(xy[nvtr[i][2]][0]-xy[nvtr[i][1]][0])/w;

alfa_mat[1][0]=(xy[nvtr[i][2]][0]*xy[nvtr[i][0]][1]-xy[nvtr[i][0]][0]*xy[nvtr[i][2]][1])/w;

alfa_mat[1][1]=(xy[nvtr[i][2]][1]-xy[nvtr[i][0]][1])/w;

alfa_mat[1][2]=(xy[nvtr[i][0]][0]-xy[nvtr[i][2]][0])/w;

alfa_mat[2][0]=(xy[nvtr[i][0]][0]*xy[nvtr[i][1]][1]-xy[nvtr[i][1]][0]*xy[nvtr[i][0]][1])/w;

alfa_mat[2][1]=(xy[nvtr[i][0]][1]-xy[nvtr[i][1]][1])/w;

alfa_mat[2][2]=(xy[nvtr[i][1]][0]-xy[nvtr[i][0]][0])/w;

 

}

//матрица жёсткости

void local_matr_gest(int num_el)

{

int q=lyam();

double w=fabs(det_D(num_el));

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

{

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

{

local_gest[i][j]=q*(w/2.0)*(alfa_mat[i][1]*alfa_mat[j][1]+alfa_mat[i][2]*alfa_mat[j][2]);

}

}

}

 

//матрица массы

void local_matr_mass(int number_obl, int num_el)

{

int m1[3][3], m2[3][3], m3[3][3];

double w=fabs(det_D(num_el));

int q=gamma(number_obl);

m1[0][0]=m2[1][1]=m3[2][2]=6.0;

m1[0][1]=m1[0][2]=m1[1][0]=m1[1][1]=m1[2][0]=m1[2][2]=m2[0][0]=m2[0][1]=m2[1][0]=m2[1][2]=m2[2][1]=m2[2][2]=m3[0][0]=m3[0][2]=m3[1][1]=m3[1][2]=m3[2][0]=m3[2][1]=2.0;

m1[1][2]=m1[2][1]=m2[0][2]=m2[2][0]=m3[0][1]=m3[1][0]=1.0;

 

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

{

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

{

local_mas[i][j]=(q*w/120.0)*(m1[i][j]+m2[i][j]+m3[i][j]);         

}

}

}

 

//локальный вектор правых частей

void lokal_vector(int number_obl, int num_el)

{

double b1[3], b2[3], b3[3];

b1[0]=b2[1]=b3[2]=2.0;

b1[1]=b1[2]=b2[0]=b2[2]=b3[0]=b3[1]=1.0;

double v1=koef_f(xy[nvtr[num_el][0]][0],xy[nvtr[num_el][0]][1]);

double v2=koef_f(xy[nvtr[num_el][1]][0],xy[nvtr[num_el][1]][1]);

double v3=koef_f(xy[nvtr[num_el][2]][0],xy[nvtr[num_el][2]][1]);

 

double w=fabs(det_D(num_el));

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

{

local_f[i]=(w/24.0)*(v1*b1[i]+v2*b2[i]+v3*b3[i]);

}

}

 

//локальная матрица А

void local_matrica()

{

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

{

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

{

local[i][j]=local_gest[i][j]+local_mas[i][j];

}

}

}

 

//глобальная матрица и вектор

void global_all(int number_el)

{

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

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

{

if(nvtr[number_el][j]<nvtr[number_el][i])

{

int h=ig[nvtr[number_el][i]];

for(int k=0;k==0;h++)

{

if(jg[h]==nvtr[number_el][j])

{

gg[h]+=local[i][j];

k=1;

}

}

}

}

for(int i=0;i<3;i++) diag[nvtr[number_el][i]]+=local[i][i];

for(int i=0;i<3;i++) glob_f[nvtr[number_el][i]]+=local_f[i];

}

double U_g(int uzl)

{

double x=xy[uzl][0];

double y=xy[uzl][1];

double Ug=x*x;

return Ug;

}

// первые кравые условия

void per_kra()

{

FILE *f;

int k,i_1,i_2,num_el;

double Ug;

f=fopen("kr_1.txt", "r");

while(!feof(f))

{

fscanf(f,"%d", &k);//номер узла

diag[k]=C;

Ug=U_g(k);

glob_f[k]=C*Ug;

}

fclose(f);

}

// вторые краевые условия

void vtor_kra()

{

FILE *f;

int i1,i2,i;

double mes,teta[2],CT[2];

    f=fopen("kr_2.txt", "r");

while(!feof(f))

{

fscanf(f,"%d %lf", &i1,&teta[0]);// номер узла, соответствующая ему тетта

fscanf(f,"%d %lf", &i2,&teta[1]);

mes=dlina(i1,i2);

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

CT[i]=0.0;

CT[0]=(2*teta[0]+teta[1])*(mes/6.0);

CT[1]=(2*teta[1]+teta[0])*(mes/6.0);

 

glob_f[i1]+=CT[0];

glob_f[i2]+=CT[1];

}

fclose(f);

}

//третьи кравевые условия

 void tret_kra()

 {

FILE *f;

int i1,i2,i,j;

double mes,C1[2][2],C12[2][2],CT[2],beta,ub[2];

    C1[0][0]=2.0; C1[0][1]=1.0;

    C1[1][0]=1.0; C1[1][1]=2.0;

f=fopen("kr_3.txt", "r");

while(!feof(f))

{

fscanf(f,"%lf", &beta);

fscanf(f,"%d %lf", &i1,&ub[0]);

fscanf(f,"%d %lf", &i2,&ub[1]);

mes=dlina(i1,i2);

for(i=0;i<2;i++) CT[i]=0.0;

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

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

C12[i][j]=(C1[i][j]*mes/6.)*beta;  

int k=1;

for(i=0; i<i2 && k;i++)

{ 

if(jg[ig[i2]+i]==i1)

{ 

k=0;

gg[ig[i2]+i]+=C12[0][1]; //учитываем матрицу С12 в глобальной матрице

}

}

diag[i1]+=C12[0][0];

diag[i2]+=C12[1][1];

 

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

{

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

{

CT[i]+=C12[i][j]*ub[j]; //умножаем матрицу на вектор (ub1, ub2)

}

}

 

glob_f[i1]+=CT[0]; //учёт вектора CT в локальном векторе правой части

glob_f[i2]+=CT[1];

}

fclose(f);

}

 

 

//головная программа

void main()

{

FILE *f3,*f4,*f5,*f6;

f3=fopen("ig.txt","r");

f4=fopen("jg.txt","r");

f5=fopen("diag.txt","w+");

f6=fopen("gg.txt","w+");

 

vvod();//считываем данные

vvod_koord();

glob_f =new double [kol_uzl];

ig = new int [kol_uzl+1]; 

diag = new double [kol_uzl];

x=new double[kol_uzl];

r=new double[kol_uzl];

p=new double[kol_uzl];

z=new double[kol_uzl];

q=new double[kol_uzl];

diag_=new double[kol_uzl];

f=new double[kol_uzl];

analitika(f);

 

for(int i=0;i<kol_uzl+1;i++)

{

fscanf(f3,"%d", &ig[i]);

}

int k=ig[kol_uzl];

 

jg =new int [k];

gg = new double [k];

gg_=new double[k];

 

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

{

fscanf(f4,"%d", &jg[i]);

}

 

null_all();

 

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

{

int number_obl=nvkat[i];

alfa(i);

local_matr_gest(i);

local_matr_mass(number_obl,i);

local_matrica();

lokal_vector(number_obl,i);

global_all(i);

}

tret_kra();//учет краевых условий

vtor_kra();

per_kra();

 

for(int i=0;i<kol_uzl;i++) fprintf(f5,"%.10ef\n", diag[i]);

for(int i=0;i<k;i++) fprintf(f6,"%.10ef\n", gg[i]);

l_mat();

los_hol(x,glob_f);

fclose(f3);

fclose(f4);

delete(x);

delete(r);

delete(p);

delete(z);

delete(q);

delete(f);

delete(diag);

delete(diag_);

delete(ig);

delete(jg);

delete(gg);

delete(gg_);

delete(glob_f);

fclose(f5);

fclose(f6);

}

 

 

 

 

Исследования :

Ω=[2;4]*[5;9]

4 конечных элемента, 5 узлов

Тест 1:     

Аналитическое решение