Прога для 425гр (мм)

erpel



#include <conio.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
typedef struct vect{
double* data;
int size;
} str_vect;
typedef str_vect* Vect;
typedef struct matrix{
double* data;
int size;
} str_matrix;
typedef str_matrix* Matrix;
int nf=1;
int N=2*2*2;
int M=(N+1)*(N+2)/2;
double h=1./N;
double d0=10,d1=1,d2=1,g1=1,g2=1,g3=1;
//double (*pFuncint,int);
void ProductA(Vect ,Vect result);
void Alg_5(Vect b,Vect* pc);
double Func(int number,int a,int b);
double Func(int number,int a,int b)
{
double x,y;
x = -0.5 + b*h +0.5*a*h;
y = (sqrt(3.)/2)*a*h;
switch(number)
{
case 1:return d0;
case 2:return d1*x+d2*y;
case 3:return g1*x*x+g2*y*y+g3*x*y;
}
return 0;
}
Vect CreateVect(int n)
{
int i;
Vect v;
v = (Vect)malloc(sizeof(str_vect;
v->size = n;
v->data = (double*)malloc(n*sizeof(double;
for(i=0;i<n;i++)
v->data[i] = 0;
return v;
}
Vect CreateVectData(int n,double* d)
{
Vect v;
v = (Vect)malloc(sizeof(str_vect;
v->size = n;
v->data = d;
return v;
}
void DeleteVect(Vect v)
{
free(v->data);
free(v);
}
void Equal(Vect v,const Vect w)
{
int i;
for(i=0;i<v->size;i++)
v->data[i]=w->data[i];
}
double Scal(const Vect v,const Vect w)
{
int i;
double d=0;
for (i=0;i<v->size;i++)
d += (v->data[i])*(w->data[i]);
return d;
}
void Plus(const Vect v,const Vect w,Vect res)
{
int i;
for (i=0;i<v->size;i++)
res->data[i] = v->data[i] + w->data[i];
return ;
}
void Minus(const Vect v,const Vect w,Vect res)
{
int i;
//Vect res;
//CreateVect(v->size,res)
for (i=0;i<v->size;i++)
res->data[i] = v->data[i] - w->data[i];
return ;
}
void PrintVect(const Vect v)
{
int i;
for(i=0;i<v->size;i++)
printf("%lf ",v->data[i]);
printf("\n");
}
void Multi(const double a,const Vect v,Vect res)
{
int i;
for(i=0;i<v->size;i++)
res->data[i]=a*v->data[i];
return ;
}
void Minus_Multi(const Vect u,const double a,const Vect v,Vect res)
{
int i;
for(i=0;i<v->size;i++)
res->data[i]=u->data[i]-a*v->data[i];
return ;
}
void Plus_Plus_Multi(Vect c,const double a1,const Vect v1,const double a2,const Vect v2)
{
int i;
for(i=0;i<c->size;i++)
c->data[i]+=a1*v1->data[i]+a2*v2->data[i];
return ;
}
int BICGTAB(const Vect b,Vect res)
{
int n,i,j ;
Vect r,r1,c,c1,p,p1,s,z,v1,v2,v3,ap,as;
double a=0,w=0,eps,g=0;
double *d;
n = b->size;
d = (double*)malloc(n*sizeof(double;
for(i=0;i<n;i++) d[i]=1;

z =CreateVectData(n,d);
r =CreateVect(n);
r1=CreateVect(n);
c =CreateVect(n);
c1=CreateVect(n);
p =CreateVect(n);
p1=CreateVect(n);
s =CreateVect(n);
v1=CreateVect(n);
v2=CreateVect(n);
v3=CreateVect(n);
ap=CreateVect(n);
as=CreateVect(n);

ProductA(c,v1);
Minus(b,v1,r);
Equal(p,r);
eps = sqrt(Scal(r,r/10000;

for(j=0; ;j++)
{
ProductA(p,ap);
a = Scal(r,z)/Scal(ap,z);

Multi(a,ap,v1);
Minus(r,v1,s);

//if (Scal(s,s)==0) {printf("Error\n");break;}
ProductA(s,as);
w = Scal(as,s)/Scal(as,as);

Multi(a,p,v1);
Plus(c,v1,v2);
Multi(w,s,v3);
Plus(v2,v3,c1);

Multi(w,as,v1);
Minus(s,v1,r1);

g = Scal(r1,z)/Scal(r,z)*(a/w);
Multi(w,ap,v1);
Minus(p,v1,v2);
Multi(g,v2,v3);
Plus(r1,v3,p1);

if (sqrt(Scal(r,r < eps) break;
Equal(r,r1);
Equal(c,c1);
Equal(p,p1);
printf("%d\n",j);
}
Equal(res,c1);
DeleteVect(z);
DeleteVect(r);
DeleteVect(r1);
DeleteVect(c);
DeleteVect(c1);
DeleteVect(p);
DeleteVect(p1);
DeleteVect(s);
DeleteVect(v1);
DeleteVect(v2);
DeleteVect(v3);
DeleteVect(ap);
DeleteVect(as);
return j;
}
void ProductA(Vect x,Vect y)
{
int a=0,b=0,i=0,j=-1,k=N+1;
//первый слой
y->data[i] =2*x->data[i] + x->data[i+1]
+ x->data[k];
for(b=1;b<N;b++)
{
i++;
k++;
y->data[i] =(1*x->data[i-1]+6*x->data[i]+1*x->data[i+1])
+(2*x->data[k-1]+2*x->data[k]);
}
i++;
k++;
y->data[i] =(1*x->data[i-1]+2*x->data[i])
+(1*x->data[k-1]);

//средние слои
for(a=1;a<N;a++)
{
i++;
j++;
y->data[i] = (1*x->data[j]+2*x->data[j+1])+
+(6*x->data[i]+2*x->data[i+1])
+(1*x->data[k]);
for(b=1;b<N-a;b++)
{
i++;
j++;
k++;
y->data[i] = (2*x->data[j]+2*x->data[j+1])+
+(2*x->data[i-1]+12*x->data[i]+2*x->data[i+1])
+(2*x->data[k-1]+2*x->data[k]);

}
i++;
j++;
k++;
y->data[i] = (2*x->data[j]+1*x->data[j+1])+
+(2*x->data[i-1]+6*x->data[i])
+(1*x->data[k-1]);
j++;
}
//последний слой
i++;
j++;
y->data[i] = (1*x->data[j]+1*x->data[j+1])+
+(2*x->data[i]);
return ;
}
Vect CreateVectB(int size,int number)
{
int a,b,i=0;
Vect B;
double *data;
B=(str_vect*)malloc(sizeof(str_vect;
data = (double*)malloc(size*sizeof(double;
B->data=data;
B->size=size;

data[0] = 1*h*Func(number,0,0)/sqrt(sqrt(3;
for(b=1;b<N;b++)
{
i++;
data[i] = 3*h*Func(number,0,b)/sqrt(sqrt(3;
}
i++;
data[i] = 1*h*Func(number,0,N)/sqrt(sqrt(3;
for(a=1;a<N;a++)
{
i++;
data[i] = 3*h*Func(number,a,0)/sqrt(sqrt(3;
for(b=1;b<N-a;b++)
{
i++;
data[i] = 6*h*Func(number,a,b)/sqrt(sqrt(3;
}
i++;
data[i] = 3*h*Func(number,a,N-a)/sqrt(sqrt(3;
}
i++;
data[i] = 1*h*Func(number,N,0)/sqrt(sqrt(3;
if (size!=i+1) printf("Error");
return B;
}
double Mistake(Vect b,Vect c,int number)
{
Vect ac;
ac=CreateVect(M);
ProductA(c,ac);
printf("acc=%lf\n",Scal(ac,c;
printf("bc=%lf\n",Scal(b,c;
double ff;
switch(number)
{
case 1:
ff=d0*d0*(sqrt(3)/4.);
break;
case 2:
ff=(d1*d1 + d2*d2*3)/(32*sqrt(3;
break;
case 3:
ff=(2*g1*g1+2*g1*g2+g3*g3+18*g2*g2)*sqrt(3)/(1920);
break;
}
printf("ff=%lf\n",ff);
return sqrt(fabs(ff-2*Scal(b,c)+Scal(ac,c;
}
void main
{

Vect c,b;
c = CreateVect(M);
b = CreateVectB(M,nf);
printf("Iteration=%d\n",BICGTAB(b,c;
printf("Mistake=%lf\n",Mistake(b,c,nf;
DeleteVect(c);
DeleteVect(b);
}
/*void Alg_5(Vect b,Vect* pc)
{
int i=0,iter=0;
double *d;
Vect r,s,p,z,ap,c,as,v1,v2,v3;

c=CreateVect(M);
r=CreateVect(M);
p=CreateVect(M);
s=CreateVect(M);
ap=CreateVect(M);
as=CreateVect(M);
v1=CreateVect(M);
v2=CreateVect(M);
v3=CreateVect(M);
d = (double*)malloc(M*sizeof(double;
//for(i=0;i<M;i++) d[i]=1;
z =CreateVect(M);
z->data[0]=1;

double alfa=0.,w=0.,betta=0.;
double p0,r0,stop = 1000,scalrz;

Equal(r,b);
Equal(p,b);
r0=sqrt(Scal(r,r;
p0=r0;
do
{
ProductA(p,ap);
alfa=Scal(r,z)/Scal(ap,z);
Minus_Multi(r,alfa,ap,s);

ProductA(s,as);
w=Scal(as,s)/Scal(as,as);
scalrz = Scal(r,z); //save S(r,z)
Minus_Multi(s,w,as,r);
Plus_Plus_Multi(c,alfa,p,w,s);

betta = (Scal(r,z)/scalrz)*(alfa/w);
Minus_Multi(p,w,ap,v1);
Plus_Multi(r,betta,v1,p);

r0=sqrt(Scal(r,r;
printf("mistake =%e ,i=%d\n" ,r0/p0,iter);

iter++;
} while ( (p0/r0)<stop );
*pc=c;
return;
}*/

Majesty_2_Info

в чем прикол?

nasteniw

плиз, запость то же самое в тегах код со всем форматированием
и поясни, в чем маза. что за задача...
может, у тебя и другие есть?

Neonka

да это я для себя от прогеров постил
задание: "Наилучшее приближение функций в L2", область - равносторонний треугольник и т.д. вообщем =)
Оставить комментарий
Имя или ник:
Комментарий: