СӨЖ-3 Екіөлшемді толқын теңдеуі үшін кері есепті түйіндес градиенттер әдісімен шешіңіз: (1) (2) (4) (5) (3) Қосымша ақпарат: (6), m – тізімдегі реттік номер. Есеп беруде: 1)түйіндес градиенттер әдісімен шешу алгоритмі



бет2/6
Дата04.11.2019
өлшемі0,87 Mb.
#51156
1   2   3   4   5   6
Байланысты:
Жаңаш Шолпан СӨЖ 7 апта


Түйіндес градиенттер әдісі (14)

Беріледі

  1. Меншіктейміз

  2. 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




Достарыңызбен бөлісу:
1   2   3   4   5   6




©engime.org 2024
әкімшілігінің қараңыз

    Басты бет