Private Sub Command1_Click()
Dim ks As Single
Dim a As Single
Dim beta As Single
Dim r1 As Single
Dim n As Single
Dim cx(1 To 10, 1 To 100) As Single
Dim lx(1 To 10, 1 To 100) As Single
Dim lp(1 To 10, 1 To 100) As Single
Dim cp(1 To 10, 1 To 100) As Single
Dim k(1 To 100) As Single
Dim h As Single
Dim k1 As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim c1 As Single
Dim e As Single
Dim ls As Single
Dim th(1 To 10) As Single
Dim m As Single
For m = 1 To 10
th(m) = 1 + 0.01 * m
Next
ls = 0.5
beta = 0.95
a = 0.33
r1 = 1 / beta - 1
ks = (r1 / a) ^ (1 / (a - 1))
ks = ks * ls
Debug.Print ks
h = 2 * ks / 100
For n = 1 To 100
k(n) = n * h
For m = 1 To 10
lx(m, n) = ls
cx(m, n) = th(m) * k(n) ^ a * lx(m, n) ^ (1 - a)
Next
Next
t1 = 0
Do Until t1 > 100
For n = 10 To 90
For m = 1 To 10
k1 = k(n) + th(m) * k(n) ^ a * lx(m, n) ^ (1 - a) - cx(m, n)
n1 = k1 / h
n2 = Int(n1)
n3 = n2 + 1
uc = 0
For m1 = 1 To 10
c1 = cx(m1, n2) + (n1 - n2) * (cx(m1, n3) - cx(m1, n2))
l1 = lx(m1, n2) + (n1 - n2) * (lx(m1, n3) - lx(m1, n2))
r1 = th(m1) * a * k1 ^ (a - 1) * l1 ^ (1 - a)
uc = uc + (beta * (1 + r1)) / c1
Next
uc = 0.1 * uc
cp(m, n) = 1 / uc
w1 = th(m) * (1 - a) * k(n) ^ a * lx(m, n) ^ (-a)
lp(m, n) = 1 - cx(m, n) / w1
Next
Next
e = 0
For m = 1 To 10
For n = 10 To 90
e = e + (cx(m, n) - cp(m, n)) ^ 2 + (lx(m, n) - lp(m, n)) ^ 2
Next
Next
For m = 1 To 10
For n = 10 To 90
lx(m, n) = lp(m, n)
cx(m, n) = cp(m, n)
Next
Next
If e < 10 ^ (-5) Then t1 = 1000
t1 = t1 + 1
Debug.Print t1, e
Loop
Dim px(1 To 10, 1 To 100) As Single
Dim ps(1 To 10, 1 To 100) As Single
Dim pi As Single
Dim ix As Single
Dim p1 As Single
Dim p2 As Single
Dim p3 As Single
Dim ms As Single
ms = 10
For m = 1 To 10
For n = 1 To 100
px(m, n) = 1
Next
Next
t3 = 0
Do Until t3 > 1000
For m = 1 To 10
For n = 10 To 90
p1 = 0.9 * px(m, n)
p2 = 1.1 * px(m, n)
k1 = k(n) + th(m) * k(n) ^ a * lx(m, n) ^ (1 - a) - cx(m, n)
n1 = k1 / h
n2 = Int(n1)
n3 = n2 + 1
z1 = 0
For m1 = 1 To 10
pz = px(m1, n2) + (n1 - n2) * (px(m1, n3) - px(m1, n2))
pi = pz / p1 - 1
l1 = lx(m1, n2) + (n1 - n2) * (lx(m1, n3) - lx(m1, n2))
r1 = th(m1) * a * k1 ^ (a - 1) * l1 ^ (1 - a)
ix = (1 + r1) * (1 + pi) - 1
z1 = z1 + ms * ix / (cx(m, n) * (1 + ix)) - p1
Next
z1 = 0.1 * z1
t2 = 0
Do Until t2 > 100
z2 = 0
For m1 = 1 To 10
pz = px(m1, n2) + (n1 - n2) * (px(m1, n3) - px(m1, n2))
pi = pz / p2 - 1
l1 = lx(m1, n2) + (n1 - n2) * (lx(m1, n3) - lx(m1, n2))
r1 = th(m1) * a * k1 ^ (a - 1) * l1 ^ (1 - a)
ix = (1 + r1) * (1 + pi) - 1
z2 = z2 + ms * ix / (cx(m, n) * (1 + ix)) - p2
Next
z2 = 0.1 * z2
p3 = p2 - z2 * (p2 - p1) / (z2 - z1)
p1 = p2
p2 = p3
z1 = z2
If z2 ^ 2 < 10 ^ (-5) Then t2 = 1000
Loop
ps(m, n) = p2
Next
Next
e = 0
For m = 1 To 10
For n = 10 To 90
e = e + (px(m, n) - ps(m, n)) ^ 2
Next
Next
For m = 1 To 10
For n = 10 To 90
px(m, n) = ps(m, n)
Next
Next
If e < 10 ^ (-5) Then t3 = 10000
t3 = t3 + 1
Debug.Print t3, e
Loop
End Sub
最終更新:2009年04月24日 23:21