pro0904 - プログラム倉庫

#include <stdio.h>
#include <math.h>
int main(void)
{
double a,beta,ks,h,k1,r1,c1,e,n1,ms;
double k[101],cx[101],cp[101],px[101],ps[101];
double p1,p2,p3,i1,pi,px1;
double z1,z2;
int n,n2,n3,t,t2;
a=0.33;
beta=0.95;
ks=pow((1/beta-1)/a,1/(a-1));
h=2*ks/100;
for (n=1;n<=100;n++){
k[n] = n*h;
cx[n]=pow(k[n],a);
}
t=0;
while (t<100){
for (n=10;n<=90;n++){
k1=k[n]+pow(k[n],a)-cx[n];
n1=k1/h;
n2=floor(n1);
n3=n2+1;
c1=cx[n2]+(n1-n2)*(cx[n3]-cx[n2]);
r1=a*pow(k1,a-1);
cp[n]=c1/(beta*(1+r1));
}
e=0;
for (n=10;n<=90;n++){
e=e+pow(cx[n]-cp[n],2);
}
for (n=10;n<=90;n++){
cx[n]=cp[n];
}
if (e < 0.0001) t=1000;
t=t+1;
}
ms=20;
for (n=1;n<=100;n++){
px[n]=1;
}
t2=0;
while (t2<100){
for (n=10;n<=90;n++){
p1=0.95*px[n];
p2=1.05*px[n];
k1=k[n]+pow(k[n],a)-cx[n];
n1=k1/h;
n2=floor(n1);
n3=n2+1;
c1=cx[n2]+(n1-n2)*(cx[n3]-cx[n2]);
r1=a*pow(k1,a-1);
px1=px[n2]+(n1-n2)*(px[n3]-px[n2]);
pi=px1/p1-1;
i1=(1+pi)*(1+r1)-1;
z1=beta*ms*i1/(c1*(1+pi))-p1;
t=0;
while (t<100){
pi=px1/p2-1;
i1=(1+pi)*(1+r1)-1;
z2=beta*ms*i1/(c1*(1+pi))-p2;
p3=p2-z2*(p2-p1)/(z2-z1);
p1=p2;
p2=p3;
z1=z2;
if (z2*z2<0.0001) t=1000;
t=t+1;
}
ps[n]=p2;
}
e=0;
for (n=10;n<=90;n++){
e=e+pow(px[n]-ps[n],2);
}  
for (n=10;n<=90;n++){
px[n]=ps[n];
}
if (e<0.0001) t2=1000;
t2=t2+1;
}
printf("%f",e);
return 0;
}