Так как вместо решения бесконечномерной задачи ищется приближённое решение, полученное как линейная комбинация финитных кусочно-полиномиальных функций, то мы получаем приближённое решение. В данной работе искомая функция и правая часть представлены в виде билинейных интерполянтов, соответственно полученные решения вычисляются точно для полиномов первой степени, что подтверждается результатами тестов. Для полиномов выше первой степени решение вычисляется уже с погрешностью.
Тексты основных модулей программ
Листинг программы
#include "solver.h" //решение СЛАУ
#include "paint.h" //визуализация решения
#include "task.h" //постановка задач
//глобальная матрицы в разреженном строчном формате и вектор
int *ig, *jg, *ia, kk;
double *di, *ggl, *ggu, *gB;
double *gr, *gphi, *hr, *hphi, *q; //сетка
int **xyk, **coord, **cnd1, **cnd2, **cnd3;
int nr, nphi, L, N, K, size, n1, n2, n3;
void AddElemenit(int i, int j, double a);
void nforming(int);
double numIn(double r1, double phi1, double r2, double phi2, int i, int j, int k);
//считывание/формирование сетки
void grid()
{
FILE *fi;
fopen_s(&fi, "g.txt", "r");
fscanf_s(fi, "%d", &nr);
gr = new double[nr];
for(int i = 0; i < nr; ++i)
fscanf_s(fi, "%lf", &gr[i]);
fscanf_s(fi, "%d", &nphi);
gphi = new double[nphi];
for(int i = 0; i < nphi; ++i)
fscanf_s(fi, "%lf", &gphi[i]);
fscanf_s(fi, "%d", &L);
hr = new double[L];
hphi = new double[L];
xyk = new int*[5];
for(int i = 0; i < L; ++i)
xyk[i] = new int[L];
for(int i = 0; i < L; ++i)
{
fscanf_s(fi, "%d %d %d %d %d", &xyk[i][0], &xyk[i][1], &xyk[i][2], &xyk[i][3], &xyk[i][4]);
--xyk[i][0]; --xyk[i][1]; --xyk[i][2]; --xyk[i][3]; --xyk[i][4];
}
fscanf_s(fi, "%d", &N);
coord = new int*[N];
for(int i = 0; i < N; ++i)
{
coord[i] = new int[2];
fscanf_s(fi, "%d %d", &coord[i][0], &coord[i][1]);
--coord[i][0];
--coord[i][1];
}
for(int i = 0; i < L; ++i)
{
hr[i] = gr[xyk[i][2]]-gr[xyk[i][1]];
hphi[i] = gphi[xyk[i][4]]-gphi[xyk[i][3]];
}
gB = new double[N];
for(int i = 0; i < N; ++i)
gB[i] = 0;
di = new double[N];
q = new double[N];
}
//получение краевых условий
void cond()
{
FILE *fi;
fopen_s(&fi, "cond.txt", "r");
fscanf_s(fi, "%d", &n1);
cnd1 = new int*[n1];
for(int i = 0; i < n1; ++i)
cnd1[i] = new int[6];
for(int i = 0; i < n1; ++i)
{
fscanf_s(fi, "%d %d %d %d %d %d", &cnd1[i][0], &cnd1[i][1], &cnd1[i][2], &cnd1[i][3], &cnd1[i][4], &cnd1[i][5]);
--cnd1[i][0]; --cnd1[i][1]; --cnd1[i][2]; --cnd1[i][3]; --cnd1[i][4]; --cnd1[i][5];
}
fscanf_s(fi, "%d", &n2);
cnd2 = new int*[n2];
for(int i = 0; i < n2; ++i)
cnd2[i] = new int[6];
for(int i = 0; i < n2; ++i)
{
fscanf_s(fi, "%d %d %d %d %d %d", &cnd2[i][0], &cnd2[i][1], &cnd2[i][2], &cnd2[i][3], &cnd2[i][4], &cnd2[i][5]);
--cnd2[i][0]; --cnd2[i][1]; --cnd2[i][2]; --cnd2[i][3]; --cnd2[i][4]; --cnd2[i][5];
}
fscanf_s(fi, "%d", &n3);
cnd3 = new int*[n3];
for(int i = 0; i < n3; ++i)
cnd3[i] = new int[6];
for(int i = 0; i < n3; ++i)
{
fscanf_s(fi, "%d %d %d %d %d %d", &cnd3[i][0], &cnd3[i][1], &cnd3[i][2], &cnd3[i][3], &cnd3[i][4], &cnd3[i][5]);
--cnd3[i][0]; --cnd3[i][1]; --cnd3[i][2]; --cnd3[i][3]; --cnd3[i][4]; --cnd3[i][5];
}
fclose(fi);
}
//поиск номера узла
int search(int r, int phi)
{
for(int i = 0; i < N; ++i)
if(coord[i][0] == r)
if(coord[i][1] == phi)
return i;
return -1;
}
//возвращает глобальный номер l-конечного элемента i-локального узла
int GlobalNum(int l, int i)
{
switch(i)
{
case 0: return search(xyk[l][1], xyk[l][3]);
case 1: return search(xyk[l][2], xyk[l][3]);
case 2: return search(xyk[l][1], xyk[l][4]);
case 3: return search(xyk[l][2], xyk[l][4]);
}
return -1;
}
//теперь генерация портрета матрицы
void GeneratePortrait()
{
int i = 0, j = 0, k = 0, ver[2][6];
int **vx;
vx = new int*[2];
vx[0] = new int[6*L];
vx[1] = new int[6*L];
ig = new int[N+1];
int *flags;//массив флагов для всех узлов
flags = new int[6*L];
for(; i < L*6; ++i)
flags[i] = 0; //вначале мы ещё не заносили узлы
for(i = 0; i < N + 1; ++i)
ig[i] = 0;
for(i = 0; i < N; ++i)
di[i] = 0;
//проходим по каждой области
for(i = 0; i < L; ++i)
{
ver[0][0] = GlobalNum(i, 0);
ver[1][0] = GlobalNum(i, 1);
ver[0][1] = GlobalNum(i, 0);
ver[1][1] = GlobalNum(i, 2);
ver[0][2] = GlobalNum(i, 1);
ver[1][2] = GlobalNum(i, 2);
ver[0][3] = GlobalNum(i, 0);
ver[1][3] = GlobalNum(i, 3);
ver[0][4] = GlobalNum(i, 1);
ver[1][4] = GlobalNum(i, 3);
ver[0][5] = GlobalNum(i, 2);
ver[1][5] = GlobalNum(i, 3);
for(int j = 0; j < 6; ++j)
{
int kk = i*6 + j;
vx[0][kk] = ver[0][j];
vx[1][kk] = ver[1][j];
}
}
for(i = 0; i < L*6; ++i)
for(j = i; j < L*6; ++j)
if(vx[0][i] == vx[0][j])
if(vx[1][i] == vx[1][j])
++flags[j];
for(k = 0; k < L*6; ++k)
if(1 == flags[k])
++ig[vx[1][k]+1];
for(i = 1; i < N + 1; ++i)
ig[i] = ig[i-1] + ig[i];
size = ig[N];
jg = new int[size];
int jj = 0;
sort(vx, flags, L);
for(i = 0; i < L*6; ++i)
{
if(2 == flags[i])
++jj;
jg[i-jj] = vx[0][i];
}
ggu = new double[size];
ggl = new double[size];
for(i = 0; i < size; ++i)
ggu[i] = ggl[i] = 0;
}
//теперь занесение элеменита в глобальную матрицу
void AddElemenit(int i, int j, double a)
{
int ind;
if(i == j)
{
di[i] += a;
return;
}
if(i < j)
{
for(ind = ig[j]; ind < ig[j+1]; ++ind)
if(jg[ind] == i)
break;
ggu[ind] += a;
}
else
{
for(ind = ig[i]; ind < ig[i+1]-1; ++ind)
if(jg[ind] == j)
break;
ggl[ind] += a;
}
}
//теперь формирование глобальной матрицы, ухх, наконец-то я до этого дошёл =D
void BuildGlobal()
{
int i;
for(i = 0; i < L; ++i)
//в каждой области формируем матрицу
nforming(i);
}
//теперь мы учимся применять вторые краевые условия
void second()
{
double lb[2], h, lM[2][2], s = 0;
for(int i = 0; i < n2; ++i)
{
kk = i;
if(cnd2[i][2] == cnd2[i][3])
//вдоль ребра по phi
for(int k = 0; k < 2; ++k)
lb[k] = numIn(gr[cnd2[i][2]], gphi[cnd2[i][4]], gr[cnd2[i][3]], gphi[cnd2[i][5]], k, k, 3);
else
//вдоль ребра по r
for(int k = 0; k < 2; ++k)
lb[k] = numIn(gr[cnd2[i][2]], gphi[cnd2[i][4]], gr[cnd2[i][3]], gphi[cnd2[i][5]], k, k, 4);
gB[search(cnd2[i][2], cnd2[i][4])] += lb[0];
gB[search(cnd2[i][3], cnd2[i][5])] += lb[1];
}
}
//теперь учитываем третьи магические условия
void third()
{
double h, ub[2], local_b[2], s = 0, s1;
for(int i = 0; i < n3; ++i)
{
kk = i;
if(cnd3[i][2] == cnd3[i][3])
for(int k = 0; k < 2; ++k)
{
for(int m = 0; m < 2; ++m)
{
s1 = beta*numIn(gr[cnd3[i][2]], gphi[cnd3[i][4]], gr[cnd3[i][3]], gphi[cnd3[i][5]], k, m, 6);
AddElemenit(search(cnd3[i][2], cnd3[i][k+4]), search(cnd3[i][3], cnd3[i][m+4]), s1);
}
local_b[k] = beta*numIn(gr[cnd3[i][2]], gphi[cnd3[i][4]], gr[cnd3[i][3]], gphi[cnd3[i][5]], k, k, 8);
}
else
for(int k = 0; k < 2; ++k)
{
for(int m = 0; m < 2; ++m)
{
s1 = beta*numIn(gr[cnd3[i][2]], gphi[cnd3[i][4]], gr[cnd3[i][3]], gphi[cnd3[i][5]], k, m, 7);
AddElemenit(search(cnd3[i][k+2], cnd3[i][4]), search(cnd3[i][m+2], cnd3[i][4]), s1);
}
local_b[k] = beta*numIn(gr[cnd3[i][2]], gphi[cnd3[i][4]], gr[cnd3[i][3]], gphi[cnd3[i][5]], k, k, 9);
}
gB[search(cnd3[i][2], cnd3[i][4])] += local_b[0];
gB[search(cnd3[i][3], cnd3[i][5])] += local_b[1];
}
}
void zeroman(int i, int r1, int phi1, int r2, int phi2)
{
int j = 0;
int z1 = search(r1, phi1);
int z2 = search(r2, phi2);
di[z1] = 1;
di[z2] = 1;
gB[z1] = lordC(i, gr[r1], gphi[phi1]);
gB[z2] = lordC(i, gr[r2], gphi[phi2]);
//обнуление верхней треугольной матрицы
for(; j < size; ++j)
if(jg[j] == z1 || jg[j] == z2)
ggu[j] = 0;
//обнуление нижней треугольной матрицы
for(j = ig[z1]; j < ig[z1+1]; ++j)
ggl[j] = 0;
for(j = ig[z2]; j < ig[z2+1]; ++j)
ggl[j] = 0;
}
//теперь учитываем ГЛАВНЫЕ КРАЕВЫЕ УСЛОВИЯ
void overlord()
{
for(int i = 0; i < n1; ++i)
zeroman(i, cnd1[i][2], cnd1[i][4], cnd1[i][3], cnd1[i][5]);
}
//подынтегральные функции (k - номер пфункции)
double locMatr(double rp, double phis, double r, double phi, double hr, double hphi, int i, int j, int k)
{
//локаные линейные базис-функции
double psi[4], dpsir[4], dpsphi[4];
double rl[2], phil[2];
rl[0] = (rp+hr-r)/hr;
rl[1] = (r-rp)/hr;
phil[0] = (phis+hphi-phi)/hphi;
phil[1] = (phi-phis)/hphi;
//собрание локальной матрицы массы
if(k == 0)
{
psi[0] = rl[0]*phil[0];
psi[1] = rl[1]*phil[0];
psi[2] = rl[0]*phil[1];
psi[3] = rl[1]*phil[1];
return psi[i]*psi[j]*r;
}
//вектора правой части
if(k == 2)
{
double local_f[4], lf;
local_f[0] = f(gr[xyk[kk][1]], gphi[xyk[kk][3]],kk);
local_f[1] = f(gr[xyk[kk][2]], gphi[xyk[kk][3]],kk);
local_f[2] = f(gr[xyk[kk][1]], gphi[xyk[kk][4]],kk);
local_f[3] = f(gr[xyk[kk][2]], gphi[xyk[kk][4]],kk);
psi[0] = rl[0]*phil[0];
psi[1] = rl[1]*phil[0];
psi[2] = rl[0]*phil[1];
psi[3] = rl[1]*phil[1];
lf = local_f[0] * psi[0] + local_f[1] * psi[1] + local_f[2] * psi[2] + local_f[3] * psi[3];
return lf*psi[i]*r;
}
//краевые условия второго типа по одной координате(phi)
if(k == 3)
{
double local_theta[2];
local_theta[0] = secondC(cnd2[kk][1], gr[cnd2[kk][2]], gphi[cnd2[kk][4]]);
local_theta[1] = secondC(cnd2[kk][1], gr[cnd2[kk][3]], gphi[cnd2[kk][5]]);
psi[0] = rl[0]*phil[0];
psi[1] = rl[0]*phil[1];
return (local_theta[0]*psi[0] + local_theta[1]*psi[1])*psi[i]*r;
}
//краевые условия второго типа по одной координате(r)
if(k == 4)
{
double local_theta[2];
local_theta[0] = secondC(cnd2[kk][1], gr[cnd2[kk][2]], gphi[cnd2[kk][4]]);
local_theta[1] = secondC(cnd2[kk][1], gr[cnd2[kk][3]], gphi[cnd2[kk][5]]);
psi[0] = rl[0]*phil[0];
psi[1] = rl[1]*phil[0];
return (local_theta[0]*psi[0] + local_theta[1]*psi[1])*psi[i]*r;
}
//краевые условия третьего типа по phi для матрицы
if(k == 6)
{
psi[0] = rl[0]*phil[0];
psi[1] = rl[0]*phil[1];
return psi[i]*psi[j]*r;
}
//краевые условия третьего типа по r для матрицы
if(k == 7)
{
psi[0] = rl[0]*phil[0];
psi[1] = rl[1]*phil[0];
return psi[i]*psi[j]*r;
}
//третьи краевые условия для вектора правой части по фи
if(k == 8)
{
double ub1, ub2;
ub1 = thirdC(cnd3[kk][1], gr[cnd3[kk][2]], gphi[cnd3[kk][4]]);
ub2 = thirdC(cnd3[kk][1], gr[cnd3[kk][3]], gphi[cnd3[kk][5]]);
psi[0] = rl[0]*phil[0];
psi[1] = rl[0]*phil[1];
return (ub1*psi[0] + ub2*psi[1])*psi[i]*r;
}
//третьи краевые условия для вектора правой части по r
if(k == 9)
{
double ub1, ub2;
ub1 = thirdC(cnd3[kk][1], gr[cnd3[kk][2]], gphi[cnd3[kk][4]]);
ub2 = thirdC(cnd3[kk][1], gr[cnd3[kk][3]], gphi[cnd3[kk][5]]);
psi[0] = rl[0]*phil[0];
psi[1] = rl[1]*phil[0];
return (ub1*psi[0] + ub2*psi[1])*psi[i]*r;
}
//раскладываем лямбду по биквадратичному базису и
//считаем элементы локальной матрицы жёсткости
double local_lambda[9], lam, rq[3], phiq[3], psix[9];
int l = k - 10;
local_lambda[0] = lambda(gr[xyk[l][1]], gphi[xyk[l][3]], 0);
local_lambda[1] = lambda(gr[xyk[l][1]]+hr/2, gphi[xyk[l][3]], 0);
local_lambda[2] = lambda(gr[xyk[l][2]], gphi[xyk[l][3]], 0);
local_lambda[3] = lambda(gr[xyk[l][1]], gphi[xyk[l][3]]+hphi/2, 0);
local_lambda[4] = lambda(gr[xyk[l][1]]+hr/2, gphi[xyk[l][3]]+hphi/2, 0);
local_lambda[5] = lambda(gr[xyk[l][2]], gphi[xyk[l][3]]+hphi/2, 0);
local_lambda[6] = lambda(gr[xyk[l][1]], gphi[xyk[l][4]], 0);
local_lambda[7] = lambda(gr[xyk[l][1]]+hr/2, gphi[xyk[l][4]], 0);
local_lambda[8] = lambda(gr[xyk[l][2]], gphi[xyk[l][4]], 0);
rq[0] = 2*(rl[1] - 0.5)*(rl[1] - 1);
rq[1] = -4*rl[1]*(rl[1] - 1);
rq[2] = 2*rl[1]*(rl[1] - 0.5);
phiq[0] = 2*(phil[1] - 0.5)*(phil[1] - 1);
phiq[1] = -4*phil[1]*(phil[1] - 1);
phiq[2] = 2*phil[1]*(phil[1] - 0.5);
psix[0] = rq[0]*phiq[0];
psix[1] = rq[1]*phiq[0];
psix[2] = rq[2]*phiq[0];
psix[3] = rq[0]*phiq[1];
psix[4] = rq[1]*phiq[1];
psix[5] = rq[2]*phiq[1];
psix[6] = rq[0]*phiq[2];
psix[7] = rq[1]*phiq[2];
psix[8] = rq[2]*phiq[2];
lam = 0;
for(int z = 0; z < 9; ++z)
lam += local_lambda[z]*psix[z];
//матрица жёсткости
dpsir[0] = -1/hr*phil[0];
dpsir[1] = 1/hr*phil[0];
dpsir[2] = -1/hr*phil[1];
dpsir[3] = phil[1]/hr;
dpsphi[0] = -rl[0]/hphi;
dpsphi[1] = -rl[1]/hphi;
dpsphi[2] = rl[0]/hphi;
dpsphi[3] = rl[1]/hphi;
return lam*r*(dpsir[i]*dpsir[j] + dpsphi[i]*dpsphi[j]/(r*r));
}
//численное интегрирование(кубатурная формула Гаусса)
double numIn(double r1, double phi1, double r2, double phi2, int i, int j, int k)
{
double sum = 0.0;
double hr = r2 - r1, hphi = phi2 - phi1;
double xh = r1 + r2, yh = phi1 + phi2;
double rp = (xh/2.0)+hr/(2*sqrt(3.0));
double xm = (xh/2.0)-hr/(2*sqrt(3.0));
double yp = (yh/2.0)+hphi/(2*sqrt(3.0));
double ym = (yh/2.0)-hphi/(2*sqrt(3.0));
//краевые условия второго типа по одной координате(phi)
if(k == 3)
{
sum = hphi/2*(locMatr(r1, phi1, rp, yp, 1, hphi, i, j, k) + locMatr(r1, phi1, xm, ym, 1, hphi, i, j, k));
return sum;
}
//краевые условия второго типа по одной координате(r)
if(k == 4)
{
sum = hr/2*(locMatr(r1, phi1, rp, yp, hr, 1, i, j, k) + locMatr(r1, phi1, xm, ym, hr, 1, i, j, k));
return sum;
}
//краевые условия третьего типа по phi
if(k == 6 || k == 8)
{
sum = hphi/2*(locMatr(r1, phi1, rp, yp, 1, hphi, i, j, k) + locMatr(r1, phi1, xm, ym, 1, hphi, i, j, k));
return sum;
}
//краевые условия третьего типа по r
if(k == 7 || k == 9)
{
sum = hr/2*(locMatr(r1, phi1, rp, yp, hr, 1, i, j, k) + locMatr(r1, phi1, xm, ym, hr, 1, i, j, k));
return sum;
}
sum = hr*hphi/4*(locMatr(r1, phi1, rp, yp, hr, hphi, i, j, k) + locMatr(r1, phi1, xm, yp, hr, hphi, i, j, k) +
locMatr(r1, phi1, xm, ym, hr, hphi, i, j, k) + locMatr(r1, phi1, rp, ym, hr, hphi, i, j, k));
return sum;
}
//теперь формирование локальных матриц численно
void nforming(int l)
{
double lA[4][4], lB[4], G[4][4], M[4][4];
kk = l;
for(int i = 0; i < 4; ++i)
{
for(int j = 0; j < 4; ++j)
{
M[i][j] = numIn(gr[xyk[l][1]], gphi[xyk[l][3]], gr[xyk[l][2]], gphi[xyk[l][4]], i, j, 0);
M[i][j] *= gamma(gr[xyk[l][1]], gphi[xyk[l][3]], l);
G[i][j] = numIn(gr[xyk[l][1]], gphi[xyk[l][3]], gr[xyk[l][2]], gphi[xyk[l][4]], i, j, 10+l);
lA[i][j] = G[i][j] + M[i][j];
AddElemenit(GlobalNum(l, i), GlobalNum(l, j), lA[i][j]);
}
lB[i] = numIn(gr[xyk[l][1]], gphi[xyk[l][3]], gr[xyk[l][2]], gphi[xyk[l][4]], i, i, 2);
gB[GlobalNum(l, i)] += lB[i];
}
}
void main(int argc, char **argv)
{
//считываем сетку
grid();
//считываем краевые условия всех типов
cond();
//формируем портрет матрицы
GeneratePortrait();
//собираем глобальные матрицу и вектор
BuildGlobal();
//учитываем вторые краевые условия
second();
//третьи
third();
//главные(естественные)
overlord();
Matrix *Matr = new Matrix(ig, jg, ggl, ggu, di, gB, N);
//решаем(ЛОС с LU - факторизацией)
Solver(Matr, q);
//рисуем
//pain(argc, argv, q);
Matr->~Matrix();
//ничего
printf_s("\nTo do nothing...\n");
_getch();
}
