Түйіндес градиенттер әдісі (14)
Беріледі
Меншіктейміз
n=1,2,…
Алгоритм:
1)(1)-(6) тура есебін шешу
2) (7) бойынша функционалды есептеу
3) (8)-(12 )есебін шешу
4)(13) бойынша градиентті есептеу
5) (14)формуласымен есептеу
Пргораммалық коды
#include
#include
#include
#include
using namespace std;
int Nx = 30,Ny=30;
double h = 1./Nx, b = 5./100., tau = 1e-3;
int i,j;
void Dp(double q[101][101], double un[101][101]) {
double unn[101][101], u[101][101];
int iteracia = 0;
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
unn[i][j] = 0.;
for (i = 1; i < Nx; i++)
for (j = 1; j < Ny; j++)
u[i][j] = unn[i][j] + tau*q[i][j];
for (i = 0; i <= Nx; i++) {
u[i][0] = 0.;
u[i][Nx] = 0.;
u[0][i] = 0.;
u[Nx][i] = 0.;
}
do {
iteracia++;
for (i = 1; i < Nx; i++) {
for (j = 1; j < Ny; j++) {
un[i][j] = 2.*u[i][j] - unn[i][j] + tau*tau*(u[i-1][j] - 4.*u[i][j] + u[i+1][j] + u[i][j-1] + u[i][j+1])/(h*h);
}}
for (i = 0; i <= Nx; i++) {
un[0][i] = 0.;
un[Nx][i] = 0.;
un[i][0] = 0.;
un[i][Nx] = 0.;
}
for (i = 1; i < Nx; i++) {
for (j = 1; j < Ny; j++) {
unn[i][j] = u[i][j];
u[i][j] = un[i][j];
}}
} while (iteracia*tau < 0.25);
}
int main() {
ofstream fout("sum.txt");
int iteracia = 0;
double mu[101][101], q[101][101], alpha=0.01,q_exact[101][101], u_exact[101][101];
for (i = 0; i <= Nx; i++) {
for (j = 0; j <= Ny; j++) {
q[i][j] = 9./5.*exp(-pow((i*h - 0.5), 2)-pow((j*h - 0.5), 2)/(2.*b*b));
q_exact[i][j] = q[i][j];
}}
clock_t t;
t = clock();
Dp (q, u_exact);
for (i = 0; i <= Nx; i++) {
for (j = 0; j <= Ny; j++) {
fout « q_exact[i][j] « " ";
} fout « endl;
}
double f[101][101], grad_J[101][101], func_J;
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
f[i][j] = u_exact[i][j];
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
q[i][j] = 0.0;
do {
iteracia++;
Dp (q, mu);
// if (iteracia % 100 == 0) {
//
// double norma1 = 0,norma2=0;
// for (i = 0; i <= Nx; i++) {
// for (j = 0; j <= Ny; j++) {
// norma1 += pow(u_exact[i][j] - mu[i][j], 2);
// norma2 += pow(q_exact[i][j] - q[i][j], 2);
// }
// fout « norma1«" "«norma2 « endl;
// }}
func_J = 0.;
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
func_J += pow(mu[i][j] - f[i][j], 2);
func_J *= h*h;
cout«func_J«endl;
double psi[101][101];
for (i = 0; i <= Nx; i++)
for (j = 0; j <= Ny; j++)
psi[i][j] = 2.*(f[i][j] - mu[i][j]);
Dp (psi, mu);
for (i = 0; i <= Nx; i++){
for (j = 0; j <= Ny; j++) {
grad_J[i][j]= - mu[i][j];
q[i][j] += - alpha* grad_J[i][j];
}}
Dp (q, mu);
} while (func_J > 1e-3);
cout « "iteraca sany = " « iteracia « endl;
Dp (q, mu);
for (i = 0; i <= Nx; i++) {
for (j = 0; j <= Ny; j++) {
fout « q[i][j] « " ";
}
fout « endl;
}
t = clock() - t;
printf ("It took me %d clicks (%f seconds).\n",t,((double)t)/CLOCKS_PER_SEC);
cout « "Iteration: " « iteracia « "\n";}
t=0.1
t=0.25
t=0.50
t=0.75
Достарыңызбен бөлісу: |