正の足し算のはずが負になる。

tim · #1

現在粒子の相互作用を表すプログラムをつくってるいるんですが密度はq+=m*w[j]で二重ループで計算しているんですが密度が何故か負になってしまいます。　w[j]は正なので負なるはずはないんですが。。。

コード:

#include <stdlib.h>
#include<stdio.h>
#include <GL/glut.h>
#include <math.h>

const float h=20;
#define N 3
float  g=  9.8;
const float K=1000;
const float q0=1;
const float n=1;


int width = 600, height = 600;

/* 初期値 */
float y[N] = {10.,10,10};
float Vy[N];
float t =0;
float dt =0.01;
float tmax=10.;
float x[N] = {-10,0,10};
float Vx[N];
const float pi=3.1415;
const float a1=4/(pi*powf(h,8.));
const float a2=-30/(pi*powf(h,5));
const float m=1;

/*距離の計算*/
inline void distance(float s[],float k[], float dis[][N],float dis_2[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N;j++){
				dis[i][j]=sqrtf( powf(s[j]-s[i],2)+powf(k[j]-k[i],2)); 
				dis_2[i][j]= powf(s[j]-s[i],2)+powf(k[j]-k[i],2);
		}
	}
}

/*位置ベクトルの差*/
inline void residual_vector(float s[], float k[], float res_x[][N],float res_y[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
		res_x[i][j]=s[j]-s[i];
		res_y[i][j]=k[j]-k[i];
		}
	}
}
/*重み関数*/
inline void weight(float s[] ,float k[], float w[][N] ){
	int i,j;
	float dis[N][N], dis_2[N][N];
	distance(s,k,dis ,dis_2);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			
			w[i][j]=a1*powf((h*h-dis_2[i][j]),3);
		}
	}	
}
/*勾配*/
inline void grad(float s[N], float k[N], float g_x[][N] ,float g_y[][N]){
	int i,j;
	float  dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N];
	distance(x,y,dis,dis_2);
	residual_vector( s, k, res_x, res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			
				g_x[i][j]=a2*powf((h-dis[i][j]),2)*res_x[i][j]/dis[i][j];
				g_y[i][j]=a2*powf((h-dis[i][j]),2)*res_y[i][j]/dis[i][j];
			
		}
	}

}
/*発散*/
void Laplacian(float s[N], float k[N], float L[][N]){
	int i,j;
	float  dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N];
	residual_vector( s, k, res_x, res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if(dis[i][j] >= 0.0 && dis[i][j] < h){
				L[i][j]=a2*(3*powf((h*h-dis_2[i][j]),2)-4*h*h-dis_2[i][j]);
			}
			else{
				L[i][j]=0;
			}
		}
	}
}
/*密度*/
inline void density(float x[] ,float y[] ,float q[]){
	 float w[N][N];
	 weight(x,y,w);
	 int i,j;
	 for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			q[i]+=m*w[i][j];
		 }
	 }
}
/*圧力*/
inline void  presure(float s[],float k[],float p[]){
	int i;
	 float q[N];
	density(s,k,q);
	for(i=0; i<N; i++){
		p[i]=K*(q[i]-q0);
	}
}
/*運動方程式*/
 float  *Fx(float s[] , float k[],float Vs[],float Vk[]){
	float dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N],g_x[N][N],g_y[N][N],q[N],p[N],w[N][N];
	static float Fx[N];
	distance(s,k,dis,dis_2);
	weight(s,k,w);
	residual_vector(s,k,res_x,res_y);
	grad(s,k,g_x,g_y);
	density(s,k,q);
	presure(s,k,p);
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				Fx[i]+=0;
			else if(dis[i][j]<h)
				Fx[i]+=-m*((p[i]+p[j])/(2*q[i]))*g_x[i][j];
			else 
				Fx[i]+=0;
				
		}
	}
	return Fx ;
}
 float * Fy(float s[] , float k[],float Vs[],float Vk[]){
	float dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N],p[N],g_x[N][N],g_y[N][N],q[N];
	static float Fy[N];
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				Fy[i]+=-g;
			else if(dis[i][j]<h)
				Fy[i]+=0;
			else 
				Fy[i]+=0; 	
		}
	}
	return Fy ;
}


void Calculate(float s[],float k[],float Vs[],float Vk[])
{    
    float kx2[N], kx1[N], ky1[N],ky2[N],kx3[N],ky3[N],kx4[N],ky4[N];
	int i;
	
		for(i=0;i<N;i++) {
                kx1[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky1[i]=dt**(Fy(s,k,Vs,Vk)+i);
				
				s[i]=s[i]+kx1[i]/2.0;
				k[i]=k[i]+ky1[i]/2.0;
                kx2[i]=dt**(Fx(s,k,Vs,Vk)+i);
				ky2[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]+kx2[i]/2.0-kx1[i]/2.0;
				k[i]=k[i]+ky2[i]/2.0-ky1[i]/2.0;
                kx3[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky3[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]+kx3[i]/2.0-kx2[i]/2.0;
				k[i]=k[i]+ky3[i]/2.0-ky2[i]/2.0;

                kx4[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky4[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]-kx3[i]/2.0;
				k[i]=k[i]-ky3[i]/2.0;

                Vx[i]=Vx[i]+(kx1[i]+2.0*kx2[i]+2.0*kx3[i]+kx4[i])/6.0;
				Vy[i]=Vy[i]+(ky1[i]+2.0*ky2[i]+2.0*ky3[i]+ky4[i])/6.0;
				s[i] += Vs[i]*dt;
				k[i] += Vk[i]*dt;
				*(Fx(s,k,Vs,Vk)+i)=0;
				*(Fy(s,k,Vs,Vk)+i)=0;


		}
   
}

void  particle(void){ 
	int i;
   /* axis */
  glColor3f(.2, .3, .5);
 
  glEnd();

  /* draw points */
  glColor3f( .9, .2, .2 );
  glPointSize(5.);
  glBegin( GL_POINTS );
  for(i=0;i<N;i++)
    glVertex2f( x[i], y[i]);
  glEnd();

 
}

void display(void){
  glClearColor( .0, .0, .0, .0);
  glClear( GL_COLOR_BUFFER_BIT );
  particle();
  glutSwapBuffers();
}

void idle(void){
  int i;
  float g_x[N][N],g_y[N][N],dis[N][N],dis_2[N][N],q[N],w[N][N];
	grad(x,y,g_x,g_y);
	distance( x, y,  dis, dis_2);
    Calculate(x,y,Vx,Vy);
	density(x,y,q);
	weight(x,y,w);
	printf("%e %f %e %f %f\n",a2,dis_2[0][1],g_x[0][1],q[1],w[1][1]);
  for(i=0; i<N; i++){
	
 
  
	/* 衝突 */
	if(x[i]>100.){
		x[i]=100.;	
		Vx[i]=-0.9*Vx[i];
	}
	if(y[i]<-100.){
		y[i]=-100.;
		Vy[i]=-0.9*Vy[i];
		
	  
	}
	if(x[i]<-100.){
		x[i]=-100.;
	  Vx[i]=-0.9*Vx[i];
	}
	if(y[i]>100.){
		y[i]=100;
	  Vy[i]=-0.9*Vy[i];
	}
	glutPostRedisplay();
  }
}

int main(int argc, char** argv){
  /* init gult */
  glutInit( &argc, argv );
  glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGB );
  glutInitWindowPosition( 0, 0 );
  glutInitWindowSize( width, height );
  glutCreateWindow( "particle" );

  /* call-back */
  glutDisplayFunc( display );
  glutIdleFunc( idle );

  /* 2D-setup */
  glMatrixMode( GL_PROJECTION );
  glLoadIdentity();
  gluOrtho2D( -100., 100., -100., 100.);

  glutMainLoop();
  
  return 0;
}

へにっくす · #2

C言語では、初期化していない変数はどういう数値が入っているか分かりません。
変数で宣言するときは必ず0で初期化しましょう。

コード:

float w[N][N];//この書き方は初期値が不定。
float w[N][N]={0};//これで全部が0になる
memset(w, 0, sizeof(float)*N*N);//これでもOK

tim · #3

ありがとうございます。　言われた通りにしてみたのですがしかしやはり負になってしまいました。

usao · #4

line108 : float q[N];
これの初期化も行ったけどダメ，ということでしょうか？

tim · #5

はい　その通りです。

へにっくす · #6

tim さんが書きました:ありがとうございます。　言われた通りにしてみたのですがしかしやはり負になってしまいました。

本当に？
修正したソースも載せないで、そんなこと言われてもねえ。
18行目、23行目、53行目、65行目、81行目、96行目、108行目、116行目、139行目、140行目、158行目、223行目にある変数すべて初期化してもですか？修正したと言うのなら、ちゃんとそのコードを載せてください。でないとこちらとしては他に原因をさぐる気がおきません。(^^;　

オフトピック

一つ一つ問題を解決していく、という手順は大事なのですよ。1人で納得したってダメです。

てゆーか、dis_2とかそれぞれの関数内に変数宣言してるけど、どういう意図をもってこんなに変数宣言してるんだか不明です。

かずま · #7

density() を次のようにして、結果を教えてください。

コード:

inline void density(float x[], float y[], float q[]){
    float w[N][N];
    weight(x,y,w);
    int i,j;
    int f = 0;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            q[i]+=m*w[i][j];
            printf("q[%d]=%10.3e, w[%d][%d]=%10.3e\n", i, q[i], i, j, w[i][j]);
            if (q[i] < 0) f = 1;
         }
    }
    if (f) exit(1);
}

tim · #8

修正したと言うのなら、ちゃんとそのコードを載せてください

失礼しました。
修正したコードはこちらです。

コード:

#include <stdlib.h>
#include<stdio.h>
#include <GL/glut.h>
#include <math.h>

const float h=20;
#define N 3
float  g=  9.8;
const float K=100;
const float q0=0;
const float n=1;


int width = 600, height = 600;

/* 初期値 */
float y[N] = {10.,10,10};
float Vy[N]={0};
float t =0;
float dt =0.01;
float tmax=10.;
float x[N] = {-10,0,10};
float Vx[N]={0};
const float pi=3.1415;
const float a1=4/(pi*powf(h,8.));
const float a2=-30/(pi*powf(h,5));
const float m=1;

/*距離の計算*/
inline void distance(float s[],float k[], float dis[][N],float dis_2[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N;j++){
				dis[i][j]=sqrtf( powf(s[j]-s[i],2)+powf(k[j]-k[i],2)); 
				dis_2[i][j]= powf(s[j]-s[i],2)+powf(k[j]-k[i],2);
		}
	}
}

/*位置ベクトルの差*/
inline void residual_vector(float s[], float k[], float res_x[][N],float res_y[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
		res_x[i][j]=s[j]-s[i];
		res_y[i][j]=k[j]-k[i];
		}
	}
}
/*重み関数*/
inline void weight(float s[] ,float k[], float w[][N] ){
	int i,j;
	float dis[N][N], dis_2[N][N]={0};
	distance(s,k,dis ,dis_2);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			w[i][j]=a1*powf((h*h-dis_2[i][j]),3);
		}
	}	
}
/*勾配*/
inline void grad(float s[N], float k[N], float g_x[][N] ,float g_y[][N]){
	int i,j;
	float  dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N]={0};
	distance(x,y,dis,dis_2);
	residual_vector( s, k, res_x, res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			
				g_x[i][j]=a2*powf((h-dis[i][j]),2)*res_x[i][j]/dis[i][j];
				g_y[i][j]=a2*powf((h-dis[i][j]),2)*res_y[i][j]/dis[i][j];
			
		}
	}

}
/*発散*/
void Laplacian(float s[N], float k[N], float L[][N]){
	int i,j;
	float  dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N]={0};
	residual_vector( s, k, res_x, res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if(dis[i][j] >= 0.0 && dis[i][j] < h){
				L[i][j]=a2*(3*powf((h*h-dis_2[i][j]),2)-4*h*h-dis_2[i][j]);
			}
			else
				L[i][j]=0;
			}
		}
	}
/*密度*/
inline void density(float s[] ,float k[] ,float q[]){
	 float w[N][N]={0};
	 weight(s,k,w);
	 int i,j;
	 for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			q[i]+= m*w[i][j];
		 }
	 }
}
/*圧力*/
inline void  presure(float s[],float k[],float p[]){
	int i;
	float q[N]={0};
	density(s,k,q);
	for(i=0; i<N; i++){
		p[i]=K*(q[i]-q0);
	}
}
/*運動方程式*/
 float  *Fx(float s[] , float k[],float Vs[],float Vk[]){
	 float dis[N][N],res_x[N][N],res_y[N][N],g_x[N][N],g_y[N][N],q[N],p[N]={0};
	static float Fx[N];
	residual_vector(s,k,res_x,res_y);
	grad(s,k,g_x,g_y);
	density(s,k,q);
	presure(s,k,p);
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				Fx[i]+=0;
			else if(dis[i][j]<h)
				Fx[i]+=-m*((p[i]+p[j])/(2*q[i]))*g_x[i][j];
			else 
				Fx[i]+=0;
				
		}
	}
	return Fx ;
}
 float * Fy(float s[] , float k[],float Vs[],float Vk[]){
	 float dis[N][N],res_x[N][N],res_y[N][N],dis_2[N][N],p[N],g_x[N][N],g_y[N][N],q[N]={0};
	static float Fy[N];
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				Fy[i]+=-g;
			else if(dis[i][j]<h)
				Fy[i]+=0;
			else 
				Fy[i]+=0; 	
		}
	}
	return Fy ;
}


void Calculate(float s[],float k[],float Vs[],float Vk[])
{    
	float kx2[N], kx1[N], ky1[N],ky2[N],kx3[N],ky3[N],kx4[N],ky4[N]={0};
	int i;
	
		for(i=0;i<N;i++) {
                kx1[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky1[i]=dt**(Fy(s,k,Vs,Vk)+i);
				
				s[i]=s[i]+kx1[i]/2.0;
				k[i]=k[i]+ky1[i]/2.0;
                kx2[i]=dt**(Fx(s,k,Vs,Vk)+i);
				ky2[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]+kx2[i]/2.0-kx1[i]/2.0;
				k[i]=k[i]+ky2[i]/2.0-ky1[i]/2.0;
                kx3[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky3[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]+kx3[i]/2.0-kx2[i]/2.0;
				k[i]=k[i]+ky3[i]/2.0-ky2[i]/2.0;

                kx4[i]=dt**(Fx(s,k,Vs,Vk)+i);
                ky4[i]=dt**(Fy(s,k,Vs,Vk)+i);

				s[i]=s[i]-kx3[i]/2.0;
				k[i]=k[i]-ky3[i]/2.0;

                Vx[i]=Vx[i]+(kx1[i]+2.0*kx2[i]+2.0*kx3[i]+kx4[i])/6.0;
				Vy[i]=Vy[i]+(ky1[i]+2.0*ky2[i]+2.0*ky3[i]+ky4[i])/6.0;
				s[i] += Vs[i]*dt;
				k[i] += Vk[i]*dt;
				*(Fx(s,k,Vs,Vk)+i)=0;
				*(Fy(s,k,Vs,Vk)+i)=0;


		}
   
}

void  particle(void){ 
	int i;
   /* axis */
  glColor3f(.2, .3, .5);
 
  glEnd();

  /* draw points */
  glColor3f( .9, .2, .2 );
  glPointSize(5.);
  glBegin( GL_POINTS );
  for(i=0;i<N;i++)
    glVertex2f( x[i], y[i]);
  glEnd();

 
}

void display(void){
  glClearColor( .0, .0, .0, .0);
  glClear( GL_COLOR_BUFFER_BIT );
  particle();
  glutSwapBuffers();
}

void idle(void){
  int i;
  float g_x[N][N],g_y[N][N],dis[N][N],dis_2[N][N],q[N],w[N][N]={0};
	grad(x,y,g_x,g_y);
	distance( x, y,  dis, dis_2);
    Calculate(x,y,Vx,Vy);
	density(x,y,q);
	weight(x,y,w);
	printf("%e %f %e %f %f\n",w[0][1],dis_2[0][1],g_x[0][1],q[1],w[1][1]);
  for(i=0; i<N; i++){
	
 
  
	/* 衝突 */
	if(x[i]>100.){
		x[i]=100.;	
		Vx[i]=-0.9*Vx[i];
	}
	if(y[i]<-100.){
		y[i]=-100.;
		Vy[i]=-0.9*Vy[i];
		
	  
	}
	if(x[i]<-100.){
		x[i]=-100.;
	  Vx[i]=-0.9*Vx[i];
	}
	if(y[i]>100.){
		y[i]=100;
	  Vy[i]=-0.9*Vy[i];
	}
	glutPostRedisplay();
  }
}

int main(int argc, char** argv){
  /* init gult */
  glutInit( &argc, argv );
  glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGB );
  glutInitWindowPosition( 0, 0 );
  glutInitWindowSize( width, height );
  glutCreateWindow( "particle" );

  /* call-back */
  glutDisplayFunc( display );
  glutIdleFunc( idle );

  /* 2D-setup */
  glMatrixMode( GL_PROJECTION );
  glLoadIdentity();
  gluOrtho2D( -100., 100., -100., 100.);

  glutMainLoop();
  
  return 0;
}

てゆーか、dis_2とかそれぞれの関数内に変数宣言してるけど、どういう意図をもってこんなに変数宣言してるんだか不明です。

自分としましては関数の引数に配列を渡したいのですが、配列を渡す方法が引用する関数の引数として使用する方法しか知らないからです。

box · #9

tim さんが書きました:
コード:
				Fx[i]+=0;
				Fx[i]+=0;
				Fy[i]+=0;
				Fy[i]+=0; 	

このような、「実質的に何もしていない」文を見ると、
「これは本当に意図どおりのコードなんだろうか？」という疑問がわいてしまいます。
「書かなくていいじゃん」というわけです。

それから、

tim さんが書きました:
コード:
 float  *Fx(float s[] , float k[],float Vs[],float Vk[]){
    static float Fx[N];

このように、関数名と変数名が同じ状態で、

tim さんが書きました:
コード:
    return Fx ;

こうしたとき、呼び出し元には何が返るんだろうか？という（素人の浅はかな）疑問もあります。

tim · #10

かずまさんが書きました:density() を次のようにして、結果を教えてください。

コード:

inline void density(float x[], float y[], float q[]){
    float w[N][N];
    weight(x,y,w);
    int i,j;
    int f = 0;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            q[i]+=m*w[i][j];
            printf("q[%d]=%10.3e, w[%d][%d]=%10.3e\n", i, q[i], i, j, w[i][j]);
            if (q[i] < 0) f = 1;
         }
    }
    if (f) exit(1);
}

ありがとうございます。
結果はこのようになりました
q[0]=-1.074e+008, w[0][0]=3.183e-003
q[0]=-1.074e+008, w[0][1]=1.343e-003
q[0]=-1.074e+008, w[0][2]=0.000e+000
q[1]=-1.074e+008, w[1][0]=1.343e-003
q[1]=-1.074e+008, w[1][1]=3.183e-003
q[1]=-1.074e+008, w[1][2]=1.343e-003
q[2]=-1.074e+008, w[2][0]=0.000e+000
q[2]=-1.074e+008, w[2][1]=1.343e-003
q[2]=-1.074e+008, w[2][2]=3.183e-003
続行するには何かキーを押してください . . .

tim · #11

box さんが書きました:
tim さんが書きました:
コード:
				Fx[i]+=0;
				Fx[i]+=0;
				Fy[i]+=0;
				Fy[i]+=0; 	
このような、「実質的に何もしていない」文を見ると、
「これは本当に意図どおりのコードなんだろうか？」という疑問がわいてしまいます。
「書かなくていいじゃん」というわけです。

あ　アホでした　その通りです。

box さんが書きました:
それから、
tim さんが書きました:
コード:
 float  *Fx(float s[] , float k[],float Vs[],float Vk[]){
    static float Fx[N];
このように、関数名と変数名が同じ状態で、
tim さんが書きました:
コード:
    return Fx ;
こうしたとき、呼び出し元には何が返るんだろうか？という（素人の浅はかな）疑問もあります。

なるほど　変数と関数を違う名前にします。

box · #12

tim さんが書きました: 結果はこのようになりました
q[0]=-1.074e+008, w[0][0]=3.183e-003
q[0]=-1.074e+008, w[0][1]=1.343e-003
q[0]=-1.074e+008, w[0][2]=0.000e+000
q[1]=-1.074e+008, w[1][0]=1.343e-003
q[1]=-1.074e+008, w[1][1]=3.183e-003
q[1]=-1.074e+008, w[1][2]=1.343e-003
q[2]=-1.074e+008, w[2][0]=0.000e+000
q[2]=-1.074e+008, w[2][1]=1.343e-003
q[2]=-1.074e+008, w[2][2]=3.183e-003
続行するには何かキーを押してください . . .

さて、かずまさんからのご指摘によるこの結果に関する考察はいかがでしょうか。

tim · #13

box さんが書きました:
tim さんが書きました: 結果はこのようになりました
q[0]=-1.074e+008, w[0][0]=3.183e-003
q[0]=-1.074e+008, w[0][1]=1.343e-003
q[0]=-1.074e+008, w[0][2]=0.000e+000
q[1]=-1.074e+008, w[1][0]=1.343e-003
q[1]=-1.074e+008, w[1][1]=3.183e-003
q[1]=-1.074e+008, w[1][2]=1.343e-003
q[2]=-1.074e+008, w[2][0]=0.000e+000
q[2]=-1.074e+008, w[2][1]=1.343e-003
q[2]=-1.074e+008, w[2][2]=3.183e-003
続行するには何かキーを押してください . . .
さて、かずまさんからのご指摘によるこの結果に関する考察はいかがでしょうか。

自分の計算でも同じ結果になってるんですけど、何故負になるのかさっぱり？？？というのが本音です。　どこで間違ったのかさっぱりです。　

珈琲 · #14

ブレークポイントは使っていますか？
またステップインとステップアウトは使用していますか？

へにっくす · #15

コード:

114| float dis[N][N],res_x[N][N],res_y[N][N],g_x[N][N],g_y[N][N],q[N],p[N]={0};

ねえ、これで全部0で初期化してると思ってるの？
下記のようにしないと全部の変数を初期化したことにならないよ。

コード:

114| float dis[N][N]={0},res_x[N][N]={0},res_y[N][N]={0},g_x[N][N]={0},g_y[N][N]={0},q[N]={0},p[N]={0};

53行目、64行目、80行目、114行目、115行目、135行目、136行目、154行目、219行目全てに言えることです。もう一度修正して載せてください。

usao · #16

ところで初期化云々以前の根本的な話なんですが，
>w[j]は正なので負なるはずはないんですが。。。
とのことですが，

コード:

w[i][j]=a1*powf((h*h-dis_2[i][j]),3);

この式って一般に ”負にならないようには見えない” のですけど，その点から勘違いされていたりはしないのでしょうか？

コードが全体的にとても読みにくい→読みたくない　ので，ざっと見での憶測ですけど，
dis_2[][] は，(まともに動けば)粒子間距離の２乗値なのだと思うのですが
常に ( h*h(=400) >= dis_2[][] ) なる関係が成り立つような仕組みがどこかにあるんでしょうか？

tim · #17

へにっくすさんが書きました:
コード:
114| float dis[N][N],res_x[N][N],res_y[N][N],g_x[N][N],g_y[N][N],q[N],p[N]={0};
ねえ、これで全部0で初期化してると思ってるの？
下記のようにしないと全部の変数を初期化したことにならないよ。
コード:
114| float dis[N][N]={0},res_x[N][N]={0},res_y[N][N]={0},g_x[N][N]={0},g_y[N][N]={0},q[N]={0},p[N]={0};
53行目、64行目、80行目、114行目、115行目、135行目、136行目、154行目、219行目全てに言えることです。もう一度修正して載せてください。

返信が遅れました
ご指摘の通り変数の初期化がなされていない事が原因でした。　ありがとうございました。

usao さんが書きました:ところで初期化云々以前の根本的な話なんですが，
>w[j]は正なので負なるはずはないんですが。。。
とのことですが，
コード:
w[i][j]=a1*powf((h*h-dis_2[i][j]),3);
この式って一般に ”負にならないようには見えない” のですけど，その点から勘違いされていたりはしないのでしょうか？

コードが全体的にとても読みにくい→読みたくない　ので，ざっと見での憶測ですけど，
dis_2[][] は，(まともに動けば)粒子間距離の２乗値なのだと思うのですが
常に ( h*h(=400) >= dis_2[][] ) なる関係が成り立つような仕組みがどこかにあるんでしょうか？

すいません　重み関数(ｗ[N][N]の事です)は極限でδ関数になる関数でr>hで０になるようにしていたんですけど、弄ってる内にとってしまいました。　
正常に動いたプログラムを書いておきます。　回答してくれた方ありがとうございました！

コード:

#include <stdlib.h>
#include<stdio.h>
#include <GL/glut.h>
#include <math.h>


int width = 600, height = 600;
#define N 24

/* 初期値 */
float x[N] = {-3,-1, 1, 3,-3,-1, 1, 3,-3,-1, 1, 3,-3,-1, 1, 3,-3,-1, 1, 3,-3,-1, 1, 3};
float y[N] = { 2, 2, 2, 2, 0, 0, 0, 0,-2,-2,-2,-2,-4,-4,-4,-4,-6,-6,-6,-6,-8,-8,-8,-8};
float Vx[N]={0};
float Vy[N]={0};
float t =0;
float dt =0.01;
float tmax=10.;
/*定数*/
const float pi=3.1415;
const float h=10;
const float a1=4/(pi*powf(h,8.));
const float a2=-30/(pi*powf(h,5));
const float m=100;
const float  g=  9.8;
const float K=1;
const float q0=100;
const float n=-1;
const float k=0.9;

/*距離の計算*/
inline void distance(float s[],float k[], float dis[][N],float dis_2[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N;j++){
				dis[i][j]=sqrtf( powf(s[j]-s[i],2)+powf(k[j]-k[i],2)); 
				dis_2[i][j]= powf(s[j]-s[i],2)+powf(k[j]-k[i],2);
		}
	}
}

/*位置ベクトルの差*/
inline void residual_vector(float s[], float k[], float res_x[][N],float res_y[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
		res_x[i][j]=s[j]-s[i];
		res_y[i][j]=k[j]-k[i];
		}
	}
}
/*重み関数*/
inline void weight(float dis[][N],float dis_2[][N], float w[][N] ){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if(dis[i][j]<h)
			w[i][j]=a1*powf((h*h-dis_2[i][j]),3);
		}
	}	
}
/*圧力用の勾配*/
inline void grad(float dis[][N],float res_x[][N],float res_y[][N], float g_x[][N] ,float g_y[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
				if(dis[i][j]<h){
				g_x[i][j]=a2*powf((h-dis[i][j]),2)*res_x[i][j]/dis[i][j];
				g_y[i][j]=a2*powf((h-dis[i][j]),2)*res_y[i][j]/dis[i][j];
				}
			
		}
	}
}
/*密度*/
inline void density(float w[][N], float q[]){
    int i,j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            q[i]+=m*w[i][j];
			 printf("q[%d]=%10.3e, w[%d][%d]=%10.3e\n", i, q[i], i,j, w[i][j]);
         }
    }
    
}
/*圧力*/
inline void  presure(float q[],float p[]){
	int i;
	for(i=0; i<N; i++){
		
		p[i]=K*(q[i]-q0);
		
	}
}
/*運動方程式*/
 inline float  *Fx(float Vs[],float dis[][N],float dis_2[][N],float p[N],float q[N],float g_x[][N],float g_y[][N],float res_x[][N],float res_y[][N]){
	static float fx[N];
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				fx[i]+=0;
			else if(dis[i][j]<h)
				fx[i]+=-m*((p[i]+p[j])/(q[i]*q[i]+q[j]*q[j]))*g_x[i][j]+n*m*(Vs[j]-Vs[i])*(res_x[i][j]*g_x[i][j]+res_y[i][j]*g_y[i][j])/(q[i]*q[j]*dis_2[i][j]);
				
		}
	}
	return fx ;
}
 inline float * Fy(float Vk[],float dis[][N],float dis_2[][N],float p[N],float q[N],float g_x[][N],float g_y[][N],float res_x[][N],float res_y[][N]){
	static float fy[N];
	int j,i;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j )
				fy[i]+=-g;
			else if(dis[i][j]<h)
				fy[i]+=-m*((p[i]+p[j])/(q[i]*q[i]+q[j]*q[j]))*g_y[i][j]+n*m*(Vk[j]-Vk[i])*((res_x[i][j]*g_x[i][j]+res_y[i][j]*g_y[i][j]))/(q[i]*q[j]*dis_2[i][j]);
		}
	}
	return fy ;
}


void runge(float s[],float k[],float Vs[],float Vk[],float dis[][N],float dis_2[][N],float p[N],float q[N],float g_x[][N],float g_y[][N],float res_x[][N],float res_y[][N])
{    
	float kx2[N]={0}, kx1[N]={0}, ky1[N]={0},ky2[N]={0},kx3[N]={0},ky3[N]={0},kx4[N]={0},ky4[N]={0};
	int i;
	
		for(i=0;i<N;i++) {
                kx1[i]=dt**(Fx(Vs,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);
                ky1[i]=dt**(Fy(Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);
				
				s[i]=s[i]+kx1[i]/2.0;
				k[i]=k[i]+ky1[i]/2.0;
                kx2[i]=dt**(Fx(Vs,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);
				ky2[i]=dt**(Fy(Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);

				s[i]=s[i]+kx2[i]/2.0-kx1[i]/2.0;
				k[i]=k[i]+ky2[i]/2.0-ky1[i]/2.0;
                kx3[i]=dt**(Fx(Vs,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);
                ky3[i]=dt**(Fy(Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);

				s[i]=s[i]+kx3[i]/2.0-kx2[i]/2.0;
				k[i]=k[i]+ky3[i]/2.0-ky2[i]/2.0;

                kx4[i]=dt**(Fx(Vs,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);
                ky4[i]=dt**(Fy(Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i);

				s[i]=s[i]-kx3[i]/2.0;
				k[i]=k[i]-ky3[i]/2.0;

                Vx[i]=Vx[i]+(kx1[i]+2.0*kx2[i]+2.0*kx3[i]+kx4[i])/6.0;
				Vy[i]=Vy[i]+(ky1[i]+2.0*ky2[i]+2.0*ky3[i]+ky4[i])/6.0;
				s[i] += Vs[i]*dt;
				k[i] += Vk[i]*dt;
				*(Fx(Vs,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i)=0;
				*(Fy(Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y)+i)=0;


		}
   
}
void calculate(float s[],float k[],float Vs[],float Vk[]){
	float dis[N][N]={0},dis_2[N][N]={0},res_x[N][N]={0},res_y[N][N]={0},g_x[N][N]={0},g_y[N][N]={0},q[N]={0},p[N]={0},w[N][N]={0};
	distance(s,k,dis,dis_2);
	residual_vector(s,k,res_x,res_y);
	weight(dis,dis_2,w);
	grad(dis,res_x,res_y,g_x,g_y);
	density(w,q);
	presure(q,p);
	runge(s,k,Vs,Vk,dis,dis_2,p,q,g_x,g_y,res_x,res_y);
}

void  particle(void){ 
	int i;
   /* axis */
  glColor3f(.2, .3, .5);
 
  glEnd();

  /* draw points */
  glColor3f( .9, .2, .2 );
  glPointSize(5.);
  glBegin( GL_POINTS );
  for(i=0;i<N;i++)
    glVertex2f( x[i], y[i]);
  glEnd();

 
}

void display(void){
  glClearColor( .0, .0, .0, .0);
  glClear( GL_COLOR_BUFFER_BIT );
  particle();
  glutSwapBuffers();
}

void idle(void){
  int i;
    calculate(x,y,Vx,Vy);
	
  for(i=0; i<N; i++){
	
 
	/*摩擦*/
	if(y[i]<=-100)
		Vx[i]=0.9*Vx[i];
	/* 衝突 */
	if(x[i]>100.){
		Vx[i]=-k*Vx[i]*(100/x[i]);
		x[i]=100.;
	}
	if(y[i]<-100.){
		Vy[i]=k*Vy[i]*(100/y[i]);
		y[i]=-100.;
	}
	if(x[i]<-100.){
	  Vx[i]=k*Vx[i]*(100/x[i]);
	  x[i]=-100.;
	}
	if(y[i]>100.){
	  Vy[i]=k*-Vy[i]*(100/y[i]);
	  y[i]=100;
	}
	glutPostRedisplay();
  }
}

int main(int argc, char** argv){
  /* init gult */
  glutInit( &argc, argv );
  glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGB );
  glutInitWindowPosition( 0, 0 );
  glutInitWindowSize( width, height );
  glutCreateWindow( "particle" );

  /* call-back */
  glutDisplayFunc( display );
  glutIdleFunc( idle );

  /* 2D-setup */
  glMatrixMode( GL_PROJECTION );
  glLoadIdentity();
  gluOrtho2D( -100., 100., -100., 100.);

  glutMainLoop();
  
  return 0;
}

正の足し算のはずが負になる。

正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

Re: 正の足し算のはずが負になる。

		4月 2024
日	月	火	水	木	金	土
	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30