<?php
for ( $m=1; $m<=10; $m++ ){
$th[$m] = 0.95 + 0.01 * $m;
}
$beta = 0.95;
$a = 0.33;
$ls=(1-$a)/(2-$a);
$z1 = 1/ $beta - 1;
$ks = $ls*pow($z1/$a,1/($a-1));
$h = 2 * $ks / 100;
for ($n = 1; $n<=100; $n++){
$k[$n] = $n * $h;
}
for ( $m=1; $m<=10; $m++ ){
for ( $n=1; $n<=100; $n++ ){
$lx[$m][$n] = 0.5;
$cx[$m][$n] = $th[$m] * pow($k[$n], $a) * pow($lx[$m][$n] ,1 - $a) ;
}
}
$t1 = 0;
Do {
for ( $m=1; $m<=10; $m++ ){
for ( $n=10; $n<=90; $n++ ){
$uc = 0;
for ( $s=1; $s<=10; $s++ ){
$k1 = $k[$n] + $th[$m] *pow($k[$n] , $a) *pow($lx[$m][$n] ,1 - $a) - $cx[$m][$n];
$n1 = $k1 / $h ;
$n2 = floor($n1);
$n3 = $n2 + 1;
$c1 = $cx[$s][$n2] + ($n1 - $n2) * ($cx[$s][$n3] - $cx[$s][$n2]);
$l1 = $lx[$s][$n2] + ($n1 - $n2) * ($lx[$s][$n3] - $lx[$s][$n2]);
$r1 = $th[$s] * $a * pow($k1 , $a - 1) *pow( $l1 , 1 - $a);
$uc = $uc + ($beta * (1 + $r1)) / $c1;
}
$uc = $uc / 10;
$cp[$m][$n] = 1 / $uc;
$w1 = $th[$m] * (1 - $a) * pow($k[$n] , $a) * pow($lx[$m][$n] , -$a);
$lp[$m][$n]= 1 - $c1 / $w1;
}
}
$e = 0;
for ( $m=1; $m<=10; $m++ ){
for ( $n=10; $n<=90; $n++ ){
$e = $e + pow($cx[$m][$n] - $cp[$m][$n] , 2) + pow($lx[$m][$n] - $lp[$m][$n], 2);
}
}
for ( $m=1; $m<=10; $m++ ){
for ( $n=10; $n<=90; $n++ ){
$cx[$m][$n] = $cp[$m][$n];
$lx[$m][$n] = $lp[$m][$n];
}
}
If ( $e < 0.01){
$t1=1000;
}
$t1 =$t1 +1;
}
while($t1 < 100);
$t=1;
$kt[$t]=$k[40];
for ( $t=1; $t<=200; $t++ ){
$m = rand(1, 10);
$tht[$t] = $th[$m];
$n1 = $kt[$t] / $h;
$n2 = floor($n1);
$n3 = $n2 + 1 ;
$ct[$t] = $cx[$m][$n2] + ($n1 - $n2) * ($cx[$m][$n3] - $cx[$m][$n2]);
$lt[$t] = $lx[$m][$n2] + ($n1 - $n2) * ($lx[$m][$n3] - $lx[$m][$n2]);
$kt[$t+1] = $kt[$t] + $tht[$t]*pow($kt[$t] , $a) * pow($lt[$t] ,1 - $a) - $ct[$t];
}
for ( $t=100; $t<=198; $t++ ){
print($ct[$t]);
print(",");
}
print($ct[199]);
?>
最終更新:2009年11月21日 10:19