max/msp gen~ 堀之内の積分公式　二重積分を計算する。

max/msp 堀之内の積分公式をgen~で書いて二重定積分を計算します。

∫∫(x+y*y)dxdyを計算しています。
範囲[in1 in2] [1  1+0.5x]
FNFとFNP1 FNP2を書き換えることで任意の定積分を計算できます。

gen~のcodeboxの中身です。

FNF(x,y){
 fx=x+y*y;  //ここに書く
 	return fx;
	}
	
FNP1(x){
 fx=1;
   return fx;
	}

FNP2(x){
 fx=1+0.5*x;
   return fx;
	}


Data cs(11),ee(11),h(11);
m=0;n=0;i=0;j=0;k=0;
s1=0;sx=0;sy=0;s=0;t=0;
x=0;y=0;fp1=0;fp2=0;yi=0;
n=10;

m=n+1;
a=in1;b=in2;

for(i=0;i<=n;i+=1){
	t=(2*i+1)*pi/(2*(n+1));
	poke(ee,t,i);
	poke(cs,cos(t),i);
	}

for(i=0;i<=n;i+=1){
	s1=0;
	t=peek(ee,i);
	for(j=2;j<=n;j+=2){
      s1 += ( cos(j*t))/(j*j-1);
		}
	poke(h,(1-2*s1)/m,i);
	}
sx =0;
for(i=0;i<=n;i+=1){
	x = ((b-a)*peek(cs,i)+(b+a))/2;
	fp1=FNP1(x);
	fp2=FNP2(x);
	sy=0;
	for(k=0;k<=n;k+=1){
	    y=( (fp2-fp1)*peek(cs,k)+(fp2+fp1))/2;
		sy += FNF(x,y)*peek(h,k);
	  }
	  yi = (fp2 - fp1) *sy;
	  sx += yi *peek(h,i);
	}

s = (b-a)*sx;

out1=s;
​