<?php
$beta = 0.95;
$a = 0.33;
$r1 = 1 / $beta - 1;
$ks = pow($r1 / $a , 1/($a-1));
$h = 2 * $ks / 100;
$t1 = 0;
for ( $n=1; $n<=100; $n++ ){
$k[$n] = $n * $h ;
$cx[$n] = pow($k[$n] , $a);
}
$t1=0;
do {
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) {
$t1=1000;
}
$t1 = $t1 + 1;
}
while($t1 < 100);
$ms = 10;
for ( $n=10; $n<=90; $n++ ){
$px[$n] = 1;
}
$t2 = 0;
do {
for ( $n=10; $n<=90; $n++ ){
$p1 = 0.9 * $px[$n];
$p2 = 1.1 * $px[$n];
$k1 = $k[$n] + pow($k[$n] ,$a) - $cx[$n];
$n1 = $k1 / $h;
$n2 = floor($n1);
$n3 = $n2 + 1;
$pz = $px[$n2] + ($n1 - $n2) * ($px[$n3] - $px[$n2]);
$pi = $pz / $p1 - 1;
$r1 = $a * pow($k1 ,$a-1);
$ix = (1 + $r1) * (1 + $pi) - 1;
$z1 = $ms * $ix / ($cx[$n] * (1 + $ix)) - $p1;
$t3 = 0;
do{
$pi = $pz / $p2 - 1;
$r1 = $a * pow($k1 ,$a-1);
$ix = (1 + $r1) * (1 + $pi) - 1;
$z2 = $ms * $ix / ($cx[$n] * (1 + $ix)) - $p2;
$p3 = $p2 - $z2 * ($p2 - $p1) / ($z2 - $z1);
$p1 = $p2;
$p2 = $p3;
$z1 = $z2;
If (pow($z2 , 2) < 0.0001) {
$t3 = 1000;
}
$t3 = $t3 + 1;
}
while($t3 < 100);
$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;
}
while($t2 < 100);
$t=1;
$kt[$t]=$k[40];
for ( $t=1; $t<=100; $t++ ){
$n1 = $kt[$t] / $h;
$n2 = floor($n1);
$n3 = $n2 + 1 ;
$ct[$t] = $cx[$n2] + ($n1 - $n2) * ($cx[$n3] - $cx[$n2]);
$pt[$t] = $px[$n2] + ($n1 - $n2) * ($px[$n3] - $px[$n2]);
$kt[$t+1] = $kt[$t] + pow($kt[$t] ,$a) - $ct[$t];
}
for ( $t=1; $t<=99; $t++ ){
print($pt[$t]);
print(",");
}
$t=100;
print($pt[$t]);
?>
最終更新:2009年06月15日 12:19