?? adi.cpp
字號:
void tridag(double a[101],double b[101],double c[101],
double r[101],double u[101],int n)
{
int nmax,j;
double bet,gam[101];
nmax = 100;
if (b[1] == 0.0) _c_exit();
bet = b[1];
u[1] = r[1] / bet;
for (j = 2; j<= n; j++)
{
gam[j] = c[j - 1] / bet;
bet = b[j] - a[j] * gam[j];
if (bet ==0.0) _c_exit();
u[j] = (r[j] - a[j] * u[j - 1]) / bet;
}
for (j = n - 1 ; j>=1; j--)
u[j] = u[j] - gam[j + 1] * u[j + 1];
}
void adi(double a[12][12],double b[12][12],double c[12][12],
double d[12][12],double e[12][12],double f[12][12],
double g[12][12],double u[12][12],int jmax,int k,
double alpha,double beta,double eps)
{
int next1, jj, kk, k1, j, n, l, kits, maxits;
double nrr, zero, two, half,nr,ab,disc,anormg,resid,nits;
double anrom, anorm;
double aa[101], bb[101], cc[101], rr[101], uu[101];
double psi[101][101], alph[7],bet[7],r[33],s[33][7],rfact;
jj = 100;
kk = 6;
nrr = pow(2 , (kk - 1));
maxits = 100;
zero = 0.0;
two = 2.0;
half = 0.5;
if (jmax > jj ) cout<< " 'increase jj'"<<endl;
if( k > kk - 1) cout<< " 'increase kk'"<<endl;
k1 = k + 1;
nr = pow(2 , k);
alph[1] = alpha;
bet[1] = beta;
for (j = 1; j<=k; j++)
{
alph[j + 1] = sqrt(alph[j] * bet[j]);
bet[j + 1] = half * (alph[j] + bet[j]);
}
s[1][1] = sqrt(alph[k1] * bet[k1]);
for (j = 1; j<= k; j++)
{
ab = alph[k1 - j] * bet[k1 - j];
for (n = 1; n<=pow( 2 , (j - 1)); n++)
{
disc = sqrt(pow(s[n][j] , 2) - ab);
s[2 * n][j + 1] = s[n][j] + disc;
s[2 * n - 1][j + 1] = ab / s[2 * n][j + 1];
}
}
for (n = 1; n<=nr; n++)
r[n] = s[n][k1];
anormg = zero;
for (j = 2; j<=jmax - 1; j++)
{
for (l = 2; l<=jmax - 1; l++)
{
anormg = anormg + fabs(g[j][l]);
psi[j][l]=-d[j][l]*u[j][l-1]+(r[1]-e[j][l])*u[j][l];
psi[j][l]=psi[j][l] - f[j][l] * u[j][l + 1];
}
}
nits = maxits / nr;
for (kits = 1; kits<= nits; kits++)
{
for (n = 1; n<=nr; n++)
{
if (n == nr)
next1 = 1;
else
next1 = n + 1;
rfact = r[n] + r[next1];
for (l = 2; l<= jmax - 1; l++)
{
for (j = 2; j<= jmax - 1; j++)
{
aa[j - 1] = a[j][l];
bb[j - 1] = b[j][l] + r[n];
cc[j - 1] = c[j][l];
rr[j - 1] = psi[j][l] - g[j][l];
}
tridag(aa,bb,cc,rr,uu,jmax - 2);
for (j = 2; j<=jmax - 1; j++)
psi[j][l] = -psi[j][l]+two*r[n]*uu[j-1];
}
for (j = 2; j<= jmax - 1; j++)
{
for (l = 2 ;l<= jmax - 1;l++)
{
aa[l - 1] = d[j][l];
bb[l - 1] = e[j][l] + r[n];
cc[l - 1] = f[j][l];
rr[l - 1] = psi[j][l];
}
tridag(aa,bb,cc,rr,uu,jmax - 2);
for (l = 2; l<= jmax - 1; l++)
{
u[j][l] = uu[l - 1];
psi[j][l] = -psi[j][l] + rfact * uu[l - 1];
}
}
}
anorm = zero;
for (j = 2; j<= jmax - 1; j++)
{
for (l = 2; l<= jmax - 1; l++)
{
resid = a[j][l]*u[j-1][l]+(b[j][l]+e[j][l])*u[j][l];
resid = resid + c[j][l] * u[j + 1][l];
resid = resid + d[j][l] * u[j][l - 1] * u[j][l - 1];
resid = resid + f[j][l] * u[j][l + 1] + g[j][l];
anrom = anorm + fabs(resid);
}
}
if (anorm < eps * anormg) return;
}
cout<< " maxits exceeded"<<endl;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -