Commit 8c5c5b96 authored by Pedro Flores's avatar Pedro Flores
Browse files

avanços

parent 520d5c64
No preview for this file type
......@@ -9,6 +9,15 @@ int main()
{
Edo a, c;
Edp b;
SL_Pentadiag slB;
slB.B = malloc(15*sizeof(double));
slB.D = malloc(15*sizeof(double));
slB.Di = malloc(15*sizeof(double));
slB.Ds = malloc(15*sizeof(double));
slB.Sdi = malloc(15*sizeof(double));
slB.Sds = malloc(15*sizeof(double));
slB.n = 15;
double tTotal;
a.p = &pA;
a.q = &qA;
a.r = &rA;
......@@ -40,14 +49,13 @@ int main()
b.yb = 45;
b.yc = 45;
b.yd = 100;
/* gaussSeidel (&a, Y);
prnVetor(Y, a.n);
gaussSeidel (&c, Y);
prnVetor(Y, c.n); */
setEdo(&a, 10, 0, 0, 12, 0, 0, &pA, &qA, &rA);
setEdp(&b,5, 3, 0, 0, L, 0, W, 20, 45, 45, 100, &fb);
double Y[a.n];
double Yedp[b.n*b.m];
gaussSeidelEDP (&b,Yedp);
gaussSeidelEDP (&b,Yedp, &tTotal, &slB);
prnVetor(Yedp, b.n*b.m);
printf("norma L2: %f\n", normaL2Residuo(&slB, Yedp));
return 1;
}
\ No newline at end of file
No preview for this file type
......@@ -5,6 +5,63 @@
#include "utils.h"
#include "func.h"
real_t normaL2Residuo(SL_Pentadiag *SL, real_t *y)
{
double R = 0.0f;
real_t A[SL->n];
real_t res[SL->n];
for(int k = 0; k < SL->n; k++)
{
A[k] = 0.0f;
}
prnVetor(SL->B, 15);
prnVetor(y, 15);
prnVetor(SL->Sdi, 15);
prnVetor(SL->Di, 15);
prnVetor(SL->D, 15);
prnVetor(SL->Ds, 15);
prnVetor(SL->Sds, 15);
for(int i = 0; i < SL->n; i++)
{
A[i] += SL->D[i]*y[i] + SL->Di[i]*y[i] + SL->Ds[i]*y[i] + SL->Di[i]*y[i] + SL->Sdi[i]*y[i] + SL->Sds[i]*y[i];
res[i] = SL->B[i] - A[i];
}
for(int k = 0; k < SL->n; k++)
{
R = R + (res[k]*res[k]);
}
return sqrtf(R);
}
void setEdp(Edp *e, int n, int m, int id, double a, double b, double c, double d, double ya, double yb, double yc,double yd, double (* f)(double, double))
{
e->id = id;
e->f = f;
e->n = n;
e->m = m;
e->a = a;
e->b = b;
e->c = 0;
e->d = d;
e->ya = ya;
e->yb = yb;
e->yc = yc;
e->yd = yd;
}
void setEdo(Edo *e, int n, int id, double a, double b, double ya, double yb, double (* p)(double), double (* q)(double), double (* r)(double))
{
e->id = id;
e->p = p;
e->q = q;
e->r = r;
e->n = n;
e->a = a;
e->b = b;
e->ya = ya;
e->yb = yb;
}
double criterioParada(real_t *x, real_t *xPrev, int n)
{
......@@ -18,8 +75,9 @@ double criterioParada(real_t *x, real_t *xPrev, int n)
return max;
}
void gaussSeidel (Edo *edoeq, double *Y)
{
void gaussSeidel (Edo *edoeq, double *Y, double *tTotal)
{
*tTotal = timestamp();
int n = edoeq->n, i, it;
double h, xi, bi, cp, d, di, ds;
double YPrev[n];
......@@ -53,11 +111,12 @@ void gaussSeidel (Edo *edoeq, double *Y)
}
it++;
}
printf("iterações: %d\n", it);
*tTotal = timestamp() - *tTotal;
}
void gaussSeidelEDP (Edp *edpeq, double *Y)
{
/* u(i,j) são as variaveis do sistema linear não precisam ser instanciadas dentro do for*/
void gaussSeidelEDP (Edp *edpeq, double *Y, double *tTotal, SL_Pentadiag *SL)
{
*tTotal = timestamp();
int n = edpeq->n;
int m = edpeq->m;
double hx,hy, xi, yi, bi, d, di, sdi, ds, sds, erro;
......@@ -70,6 +129,7 @@ void gaussSeidelEDP (Edp *edpeq, double *Y)
for(int k = 0; k < MAXIT; k++)
{
int idx = 0;
int c = 0;
for (int j=0; j < m; ++j)
{
yi = edpeq->c + (j+1)*hy;
......@@ -79,9 +139,10 @@ void gaussSeidelEDP (Edp *edpeq, double *Y)
bi = (hx*hx)*(hy*hy) * edpeq->f(xi, yi);
sdi = hx*hx;
di = hy*hy;
d = -2 *(hx + hy);
d = -(2*hx*hx + 2*hy*hy + hx*hx*hy*hy);
ds = hy*hy;
sds = hx*hx;
// j: 0 -> m - 1 i: 0 -> n -1
// ya = 20 yb = 45; malha x
// yc = 45 yd = 100; malha y
......@@ -93,28 +154,73 @@ void gaussSeidelEDP (Edp *edpeq, double *Y)
// bi: termos indep
//di e sdi
if(j == 0 && i == 0) bi -= ds*Y[idx+1] + sds*Y[idx+2] + edpeq->yc*hy*hy + edpeq->ya*hx*hx;
if(j == 0 && i == 0){
di = 0;
sdi = 0;
bi -= ds*Y[idx+n] + sds*Y[idx+1] + edpeq->yc*hy*hy + edpeq->ya*hx*hx;
}
//di e sds
else if(j == 0 && i == n-1) bi -= ds*Y[idx+1] + sdi*Y[idx-2] + edpeq->yc*hy*hy + edpeq->yb*hx*hx;
else if(j == 0 && i == n-1){
di = 0;
sds = 0;
bi -= ds*Y[idx+n] + sdi*Y[idx-1] + edpeq->yc*hy*hy + edpeq->yb*hx*hx;
}
//ds e sdi
else if(j == m-1 && i == 0) bi -= di*Y[idx-1] + sds*Y[idx+2] + edpeq->yd*hy*hy + edpeq->ya*hx*hx;
else if(j == m-1 && i == 0) {
ds = 0;
sdi = 0;
bi -= di*Y[idx-n] + sds*Y[idx+1] + edpeq->yd*hy*hy + edpeq->ya*hx*hx;
}
//ds e sds
else if(j == m-1 && i == n-1) bi -= di*Y[idx-1] + sdi*Y[idx-2] + edpeq->yd*hy*hy + edpeq->yb*hx*hx;
else if(j == m-1 && i == n-1) {
ds = 0;
sds = 0;
bi -= di*Y[idx-n] + sdi*Y[idx-1] + edpeq->yd*hy*hy + edpeq->yb*hx*hx;
}
//di
else if (j == 0) bi -= ds*Y[idx+1] + sds*Y[idx+2] + sdi*Y[idx-2] + edpeq->yc*hy*hy;
else if (j == 0) {
// di = 0;
// bi -= ds*Y[idx+n] + sds*Y[idx+1] + sdi*Y[idx-1] + edpeq->yc*hy*hy;
bi -= (ds*Y[idx+n] + di*Y[idx-n] + sds*Y[idx+1] + edpeq->ya*hx*hx);
sdi = 0;
}
//ds
else if (j == m-1) bi -= di*Y[idx-1] + sds*Y[idx+2] + sdi*Y[idx-2] + edpeq->yd*hy*hy;
else if (j == m-1) {
bi -= (ds*Y[idx+n] + di*Y[idx-n] + sdi*Y[idx-1] + edpeq->yb*hx*hx);
sds = 0;
// bi -= di*Y[idx-n] + sds*Y[idx+1] + sdi*Y[idx-1] + edpeq->yd*hy*hy;
// ds = 0;
}
//sdi
else if(i == 0) bi -= (ds*Y[idx+1] + di*Y[idx-1] + sds*Y[idx+2] + edpeq->ya*hx*hx);
else if(i == 0) {
di = 0;
bi -= ds*Y[idx+n] + sds*Y[idx+1] + sdi*Y[idx-1] + edpeq->yc*hy*hy;
// bi -= (ds*Y[idx+n] + di*Y[idx-n] + sds*Y[idx+1] + edpeq->ya*hx*hx);
// sdi = 0;
}
//sds
else if(i == n-1) bi -= (ds*Y[idx+1] + di*Y[idx-1] + sdi*Y[idx-2] + edpeq->yb*hx*hx);
else if(i == n-1) {
// bi -= (ds*Y[idx+n] + di*Y[idx-n] + sdi*Y[idx-1] + edpeq->yb*hx*hx);
// sds = 0;
bi -= di*Y[idx-n] + sds*Y[idx+1] + sdi*Y[idx-1] + edpeq->yd*hy*hy;
ds = 0;
}
//sem ponto de contorno
else bi -= ds*Y[idx+1] + di*Y[idx-1] + sds*Y[idx+2] + sdi*Y[idx-2];
else bi -= ds*Y[idx+n] + di*Y[idx-n] + sds*Y[idx+1] + sdi*Y[idx-1];
Y[idx] = bi / d;
SL->B[idx] = bi;
SL->D[idx] = d;
SL->Di[idx] = di;
SL->Ds[idx] = ds;
SL->Sdi[idx] = sdi;
SL->Sds[idx] = sds;
idx += 1;
// prnVetor(Y, 15);
// printf("[xi: %f] [yi: %f] [j: %d, i: %d] *bi %f* *Y[%d]: %f* [c: %d]\n", xi, yi, j, i, bi, idx,Y[idx], c);
}
}
}
*tTotal = timestamp() - *tTotal;
}
// Exibe SL na saída padrão
......
// Parâmetros para teste de convergência
#define MAXIT 1 // Número máximo de iterações em métodos iterativos
#define MAXIT 50 // Número máximo de iterações em métodos iterativos
#define L 6
#define W 8
#define PI 3
typedef double real_t;
/*---------------estruturas---------------*/
// Matriz pentadiagonal
......@@ -18,6 +19,7 @@ typedef struct {
// Equação Diferencial Ordinária para MDF
typedef struct {
int n; // número de pontos internos na malha
int id; //para diferenciar quest a e c
double a, b; // intervalo
double ya, yb; // condições contorno
double (* p)(double), (* q)(double), (* r)(double);
......@@ -26,15 +28,19 @@ typedef struct {
// Equação Diferencial Parcial para MDF
typedef struct {
int n, m; // número de pontos internos na malha em x y
int id; //para diferenciar quest b e d
double a, b, c, d; // intervalo
double ya, yb, yc, yd; // condições contorno
double (* f)(double, double);
} Edp;
/*---------------funções---------------*/
void gaussSeidel (Edo *edoeq, double *Y);
void gaussSeidel (Edo *edoeq, double *Y, double *tTotal);
void prnVetor (real_t *v, unsigned int n);
void gaussSeidelEDP (Edp *edpeq, double *Y);
void gaussSeidelEDP (Edp *edpeq, double *Y, double *tTotal, SL_Pentadiag *SL);
void setEdp(Edp *e, int n, int m, int id, double a, double b, double c, double d, double ya, double yb, double yc,double yd, double (* f)(double, double));
void setEdo(Edo *e, int n, int id, double a, double b, double ya, double yb, double (* p)(double), double (* q)(double), double (* r)(double));
real_t normaL2Residuo(SL_Pentadiag *SL, real_t *y);
double pA (double num);
double qA (double num);
double rA (double num);
......
No preview for this file type
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment