# Emily Dinan, Meg Donohue, Rui Zhang
# Project 2
# Problem 1
def cobweb(g,z,x_0,n,huge):
# INPUT
# -"g" is function
# -"z" is variable defining g
# -"x_0" is initial value of z
# -"n" is number of iterations to perform
# -"huge" the value of x with absolute value bigger than this is so large that we cut short the loop
# OUTPUT
# -"P"-plot showing the iterations graphically
x = RDF(x_0)
xmin = x-RDF(0.001)
xmax = x+RDF(0.001)
P = plot(z, z, xmin, xmax, color='gray')
for i in range(n):
y = g(x);
P +=line( [ ( x, x ), ( x, y ) ], color='red')
P +=line( [ ( x, y ), ( y, y ) ], color='orange')
x = y
if (y > xmax): xmax = y
if (y < xmin): xmin = y
if (xmin < -huge): break
if (xmax > huge): break
P +=g.plot(z, xmin, xmax, color ='blue')
P +=plot(z, z, xmin, xmax, color='gray')
P.show(ymin = xmin, ymax = xmax)
return
|
|
###############
print"Project 1, A"
print"when r=2"
p=var('p')
solve([2*p*e^(-(p/1000)^2)==p],p);
Project 1, A when r=2 [p == -1000*sqrt(log(2)), p == 1000*sqrt(log(2)), p == 0] Project 1, A when r=2 [p == -1000*sqrt(log(2)), p == 1000*sqrt(log(2)), p == 0] |
p=var('p')
f(p)=diff(2*p*e^(-(p/1000)^2), p);
print (abs(f(0))-1, "unstable");print(abs(f(1000*sqrt(log(2.))))-1, "stable");
(1, 'unstable') (-0.613705638880110, 'stable') (1, 'unstable') (-0.613705638880110, 'stable') |
print"the two non-negative equilibrium are 0 and 1000*log(2)^(1/2)"
print"(1)"
cobweb(2*x*e^(-(x/1000)^2),x,700,10, 1000);
print" x=1000*log(2)^(1/2) is the stable equilibrium point";
print"(2)"
cobweb(2*x*e^(-(x/1000)^2),x,.1,10, 1000);
print" x=0 is the unstable equilibrium point";
the two non-negative equilibrium are 0 and 1000*log(2)^(1/2) (1) __main__:7: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) the two non-negative equilibrium are 0 and 1000*log(2)^(1/2) (1) __main__:7: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) |
x=1000*log(2)^(1/2) is the stable equilibrium point
(2)
__main__:12: DeprecationWarning: Substitution using function-call syntax
and unnamed arguments is deprecated and will be removed from a future
release of Sage; you can use named arguments instead, like EXPR(x=...,
y=...)
x=0 is the unstable equilibrium point############################################################
print "project 1, B"
print "when r=3"
p=var('p');
solve([3*p*e^(-(p/1000)^2)==p],p);
project 1, B when r=3 [p == -1000*sqrt(log(3)), p == 1000*sqrt(log(3)), p == 0] project 1, B when r=3 [p == -1000*sqrt(log(3)), p == 1000*sqrt(log(3)), p == 0] |
p=var('p')
f(p)=diff(3*p*e^(-(p/1000)^2), p);
print(abs(f(0))-1,"unstable");print(abs(f(1000*sqrt(log(3.))))-1,"unstable");
(2, 'unstable') (0.197224577336220, 'unstable') (2, 'unstable') (0.197224577336220, 'unstable') |
print"the two non-negative equilibrium are 0 and 1000*log(3)^(1/2)"
print"(1)"
cobweb(3*x*e^(-(x/1000)^2),x,1048,10,10000);
print" x=1000*ln(3)^(1/2) is the unstable equilibrium point";
print"(2)"
cobweb(3*x*e^(-(x/1000)^2),x,.1,10, 1000);
print" x=0 is the unstable equilibrium point";
the two non-negative equilibrium are 0 and 1000*log(3)^(1/2) (1) __main__:6: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) the two non-negative equilibrium are 0 and 1000*log(3)^(1/2) (1) __main__:6: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) |
x=1000*ln(3)^(1/2) is the unstable equilibrium point
(2)
__main__:11: DeprecationWarning: Substitution using function-call syntax
and unnamed arguments is deprecated and will be removed from a future
release of Sage; you can use named arguments instead, like EXPR(x=...,
y=...)
x=0 is the unstable equilibrium pointprint " n p"
P = vector (RDF,1000)
P[0] = 200.
for n in range(1,101): P[n]=3*P[n-1]*e^(-(P[n-1]/1000)^2);
for n in range(0,101): print "%3s %11s" %(n,P[n])
n p 0 200.0 1 576.473663491 2 1240.43829417 3 798.827626713 4 1266.01809757 5 764.668466658 6 1278.36447725 7 748.247151411 8 1282.37849113 9 742.920858546 10 1283.40298294 11 741.562525941 12 1283.64220696 13 741.245414106 14 1283.69676024 15 741.173102576 16 1283.70913141 17 741.156704522 18 1283.71193326 19 741.152990657 20 1283.71256765 21 741.152149776 22 1283.71271127 23 741.151959399 24 1283.71274379 25 741.151916297 26 1283.71275115 27 741.151906539 28 1283.71275282 29 741.15190433 30 1283.7127532 31 741.15190383 32 1283.71275328 33 741.151903717 34 1283.7127533 35 741.151903691 36 1283.71275331 37 741.151903685 38 1283.71275331 39 741.151903684 40 1283.71275331 41 741.151903684 42 1283.71275331 43 741.151903684 44 1283.71275331 45 741.151903684 46 1283.71275331 47 741.151903684 48 1283.71275331 49 741.151903684 50 1283.71275331 51 741.151903684 52 1283.71275331 53 741.151903684 54 1283.71275331 55 741.151903684 56 1283.71275331 57 741.151903684 58 1283.71275331 59 741.151903684 60 1283.71275331 61 741.151903684 62 1283.71275331 63 741.151903684 64 1283.71275331 65 741.151903684 66 1283.71275331 67 741.151903684 68 1283.71275331 69 741.151903684 70 1283.71275331 71 741.151903684 72 1283.71275331 73 741.151903684 74 1283.71275331 75 741.151903684 76 1283.71275331 77 741.151903684 78 1283.71275331 79 741.151903684 80 1283.71275331 81 741.151903684 82 1283.71275331 83 741.151903684 84 1283.71275331 85 741.151903684 86 1283.71275331 87 741.151903684 88 1283.71275331 89 741.151903684 90 1283.71275331 91 741.151903684 92 1283.71275331 93 741.151903684 94 1283.71275331 95 741.151903684 96 1283.71275331 97 741.151903684 98 1283.71275331 99 741.151903684 100 1283.71275331 n p 0 200.0 1 576.473663491 2 1240.43829417 3 798.827626713 4 1266.01809757 5 764.668466658 6 1278.36447725 7 748.247151411 8 1282.37849113 9 742.920858546 10 1283.40298294 11 741.562525941 12 1283.64220696 13 741.245414106 14 1283.69676024 15 741.173102576 16 1283.70913141 17 741.156704522 18 1283.71193326 19 741.152990657 20 1283.71256765 21 741.152149776 22 1283.71271127 23 741.151959399 24 1283.71274379 25 741.151916297 26 1283.71275115 27 741.151906539 28 1283.71275282 29 741.15190433 30 1283.7127532 31 741.15190383 32 1283.71275328 33 741.151903717 34 1283.7127533 35 741.151903691 36 1283.71275331 37 741.151903685 38 1283.71275331 39 741.151903684 40 1283.71275331 41 741.151903684 42 1283.71275331 43 741.151903684 44 1283.71275331 45 741.151903684 46 1283.71275331 47 741.151903684 48 1283.71275331 49 741.151903684 50 1283.71275331 51 741.151903684 52 1283.71275331 53 741.151903684 54 1283.71275331 55 741.151903684 56 1283.71275331 57 741.151903684 58 1283.71275331 59 741.151903684 60 1283.71275331 61 741.151903684 62 1283.71275331 63 741.151903684 64 1283.71275331 65 741.151903684 66 1283.71275331 67 741.151903684 68 1283.71275331 69 741.151903684 70 1283.71275331 71 741.151903684 72 1283.71275331 73 741.151903684 74 1283.71275331 75 741.151903684 76 1283.71275331 77 741.151903684 78 1283.71275331 79 741.151903684 80 1283.71275331 81 741.151903684 82 1283.71275331 83 741.151903684 84 1283.71275331 85 741.151903684 86 1283.71275331 87 741.151903684 88 1283.71275331 89 741.151903684 90 1283.71275331 91 741.151903684 92 1283.71275331 93 741.151903684 94 1283.71275331 95 741.151903684 96 1283.71275331 97 741.151903684 98 1283.71275331 99 741.151903684 100 1283.71275331 |
#shown by the table that after 10 trails, all P with odd n are about 741; all P with even n are about 1283.
|
|
##################
print "project 1, C"
print "when r=3.7"
p=var('p')
solve([3.7*p*e^(-(p/1000)^2)==p],p);
project 1, C when r=3.7 [p == -1000*sqrt(log(37/10)), p == 1000*sqrt(log(37/10)), p == 0] project 1, C when r=3.7 [p == -1000*sqrt(log(37/10)), p == 1000*sqrt(log(37/10)), p == 0] |
p=var('p')
f(p)=diff(3.7*p*e^(-(p/1000)^2), p);
print(abs(f(0))-1,"unstable");print(abs(f(1000*sqrt(log(37/10.))))-1,"unstable");
(2.70000000000000, 'unstable') (0.616665639300358, 'unstable') (2.70000000000000, 'unstable') (0.616665639300358, 'unstable') |
print"the two non-negative equilibrium are 0 and 1000*log(37/10)^(1/2)"
print"(1)"
cobweb(3.7*x*e^(-(x/1000)^2),x,1150,10,10000);
print" x=1000*ln(37/10)^(1/2) is the unstable equilibrium point";
print"(2)"
cobweb(3.7*x*e^(-(x/1000)^2),x,.1,10, 1000);
print" x=0 is the unstable equilibrium point";
the two non-negative equilibrium are 0 and 1000*log(37/10)^(1/2) (1) the two non-negative equilibrium are 0 and 1000*log(37/10)^(1/2) (1) |
x=1000*ln(37/10)^(1/2) is the unstable equilibrium point
(2)
x=0 is the unstable equilibrium pointprint " n p"
P = vector (RDF,1000)
P[0] = 200.
for n in range(1,101): P[n]=3.7*P[n-1]*e^(-(P[n-1]/1000)^2);
for n in range(0,101): print "%3s %11s" %(n,P[n])
n p 0 200.0 1 710.984184973 2 1586.81556012 3 473.349346683 4 1399.83115957 5 729.903586681 6 1585.23196567 7 475.258307203 8 1402.93369838 9 725.187817233 10 1585.83463125 11 474.531248222 12 1401.75511645 13 726.97769646 14 1585.62201322 15 474.787671882 16 1402.17121739 17 726.345560754 18 1585.69936161 19 474.694377343 20 1402.01988191 21 726.575440741 22 1585.67151845 23 474.727959331 24 1402.07436329 25 726.492679455 26 1585.68158007 27 474.715823708 28 1402.05467612 29 726.522585272 30 1585.67794917 31 474.720203019 32 1402.06178064 33 726.511793091 34 1585.6792601 35 474.718621877 36 1402.05921558 37 726.515689556 38 1585.67878688 39 474.71919264 40 1402.06014152 41 726.514282999 42 1585.67895771 43 474.718986591 44 1402.05980725 45 726.514790773 46 1585.67889604 47 474.719060974 48 1402.05992792 49 726.514607468 50 1585.67891831 51 474.719034122 52 1402.05988436 53 726.514673641 54 1585.67891027 55 474.719043816 56 1402.05990009 57 726.514649753 58 1585.67891317 59 474.719040316 60 1402.05989441 61 726.514658377 62 1585.67891212 63 474.71904158 64 1402.05989646 65 726.514655263 66 1585.6789125 67 474.719041124 68 1402.05989572 69 726.514656387 70 1585.67891237 71 474.719041288 72 1402.05989599 73 726.514655982 74 1585.67891241 75 474.719041229 76 1402.05989589 77 726.514656128 78 1585.6789124 79 474.71904125 80 1402.05989592 81 726.514656075 82 1585.6789124 83 474.719041242 84 1402.05989591 85 726.514656094 86 1585.6789124 87 474.719041245 88 1402.05989592 89 726.514656087 90 1585.6789124 91 474.719041244 92 1402.05989591 93 726.51465609 94 1585.6789124 95 474.719041245 96 1402.05989592 97 726.514656089 98 1585.6789124 99 474.719041245 100 1402.05989592 n p 0 200.0 1 710.984184973 2 1586.81556012 3 473.349346683 4 1399.83115957 5 729.903586681 6 1585.23196567 7 475.258307203 8 1402.93369838 9 725.187817233 10 1585.83463125 11 474.531248222 12 1401.75511645 13 726.97769646 14 1585.62201322 15 474.787671882 16 1402.17121739 17 726.345560754 18 1585.69936161 19 474.694377343 20 1402.01988191 21 726.575440741 22 1585.67151845 23 474.727959331 24 1402.07436329 25 726.492679455 26 1585.68158007 27 474.715823708 28 1402.05467612 29 726.522585272 30 1585.67794917 31 474.720203019 32 1402.06178064 33 726.511793091 34 1585.6792601 35 474.718621877 36 1402.05921558 37 726.515689556 38 1585.67878688 39 474.71919264 40 1402.06014152 41 726.514282999 42 1585.67895771 43 474.718986591 44 1402.05980725 45 726.514790773 46 1585.67889604 47 474.719060974 48 1402.05992792 49 726.514607468 50 1585.67891831 51 474.719034122 52 1402.05988436 53 726.514673641 54 1585.67891027 55 474.719043816 56 1402.05990009 57 726.514649753 58 1585.67891317 59 474.719040316 60 1402.05989441 61 726.514658377 62 1585.67891212 63 474.71904158 64 1402.05989646 65 726.514655263 66 1585.6789125 67 474.719041124 68 1402.05989572 69 726.514656387 70 1585.67891237 71 474.719041288 72 1402.05989599 73 726.514655982 74 1585.67891241 75 474.719041229 76 1402.05989589 77 726.514656128 78 1585.6789124 79 474.71904125 80 1402.05989592 81 726.514656075 82 1585.6789124 83 474.719041242 84 1402.05989591 85 726.514656094 86 1585.6789124 87 474.719041245 88 1402.05989592 89 726.514656087 90 1585.6789124 91 474.719041244 92 1402.05989591 93 726.51465609 94 1585.6789124 95 474.719041245 96 1402.05989592 97 726.514656089 98 1585.6789124 99 474.719041245 100 1402.05989592 |
#shown by the table that after 13 trails, all P[n],which n divided by 4 remains 3, are about 474; all P[n], which n can be divided by 4, are about 1402; all P[n], which n divided by 4 remains 1, are about 726; and all P[n], which n divide by 4 remains 2,are about 1585.
|
|
print"Project 1, D"
print "Using algebra to find the equilibrium points, we solved the equation"
print "x=3.8x*e^(-(x/1000)^2)and found that when r=3.8, the two equilibrium are"
print "0 and 1000*ln(3.8)^(5/19)";
Project 1, D Using algebra to find the equilibrium points, we solved the equation x=3.8x*e^(-(x/1000)^2)and found that when r=3.8, the two equilibrium are 0 and 1000*ln(3.8)^(5/19) Project 1, D Using algebra to find the equilibrium points, we solved the equation x=3.8x*e^(-(x/1000)^2)and found that when r=3.8, the two equilibrium are 0 and 1000*ln(3.8)^(5/19) |
p=var('p')
f(p)=diff(3.8*p*e^(-(p/1000)^2), p);
print(abs(f(0))-1,"unstable");print(abs(f(1000*ln(3.8)^(5/19)))-1,"unstable");
(2.80000000000000, 'unstable') (0.575853596367162, 'unstable') (2.80000000000000, 'unstable') (0.575853596367162, 'unstable') |
cobweb(3.8*x*e^(-(x/1000)^2),x,700,10, 1000);
print" x=1000*ln(3.8)^(5/19) is an unstable equilibrium point";
cobweb(3.8*x*e^(-(x/1000)^2),x,.1,10, 1000);
print" x=0 is an unstable equilibrium point";
|
x=1000*ln(3.8)^(5/19) is an unstable equilibrium point
x=0 is an unstable equilibrium pointprint " n p"
P = vector (RDF,1000)
P[0] = 200.
for n in range(1,101): P[n]=3.8*P[n-1]*e^(-(P[n-1]/1000)^2);
for n in range(0,101): print "%3s %11s" %(n,P[n])
n p 0 200.0 1 730.199973756 2 1628.03248182 3 436.871369685 4 1371.66865813 5 794.179443694 6 1606.14072191 7 462.618417235 8 1419.25565737 9 719.532916179 10 1629.25107208 11 435.466440703 12 1368.93422736 13 798.558260742 14 1603.7721042 15 465.461693929 16 1424.21531 17 711.936056349 18 1629.67553757 19 434.977769963 20 1367.97980086 21 800.088753822 22 1602.91915301 23 466.488315209 24 1425.99157197 25 709.224269171 26 1629.73678123 27 434.907292412 28 1367.84200928 29 800.309802315 30 1602.79489331 31 466.637996684 32 1426.24990726 33 708.830270398 34 1629.74170723 35 434.901624035 36 1367.83092542 37 800.32758429 38 1602.78488569 39 466.650053064 40 1426.27070824 41 708.798550287 42 1629.74205991 43 434.901218204 44 1367.83013185 45 800.328857424 46 1602.7841691 47 466.650916352 48 1426.27219764 49 708.796279081 50 1629.74208491 51 434.901189434 52 1367.8300756 53 800.328947676 54 1602.7841183 55 466.650977551 56 1426.27230322 57 708.796118075 58 1629.74208668 59 434.901187396 60 1367.83007161 61 800.32895407 62 1602.78411471 63 466.650981886 64 1426.2723107 65 708.79610667 66 1629.74208681 67 434.901187252 68 1367.83007133 69 800.328954523 70 1602.78411445 71 466.650982193 72 1426.27231123 73 708.796105862 74 1629.74208682 75 434.901187242 76 1367.83007131 77 800.328954555 78 1602.78411443 79 466.650982215 80 1426.27231127 81 708.796105804 82 1629.74208682 83 434.901187241 84 1367.83007131 85 800.328954557 86 1602.78411443 87 466.650982216 88 1426.27231127 89 708.7961058 90 1629.74208682 91 434.901187241 92 1367.83007131 93 800.328954557 94 1602.78411443 95 466.650982216 96 1426.27231127 97 708.7961058 98 1629.74208682 99 434.901187241 100 1367.83007131 n p 0 200.0 1 730.199973756 2 1628.03248182 3 436.871369685 4 1371.66865813 5 794.179443694 6 1606.14072191 7 462.618417235 8 1419.25565737 9 719.532916179 10 1629.25107208 11 435.466440703 12 1368.93422736 13 798.558260742 14 1603.7721042 15 465.461693929 16 1424.21531 17 711.936056349 18 1629.67553757 19 434.977769963 20 1367.97980086 21 800.088753822 22 1602.91915301 23 466.488315209 24 1425.99157197 25 709.224269171 26 1629.73678123 27 434.907292412 28 1367.84200928 29 800.309802315 30 1602.79489331 31 466.637996684 32 1426.24990726 33 708.830270398 34 1629.74170723 35 434.901624035 36 1367.83092542 37 800.32758429 38 1602.78488569 39 466.650053064 40 1426.27070824 41 708.798550287 42 1629.74205991 43 434.901218204 44 1367.83013185 45 800.328857424 46 1602.7841691 47 466.650916352 48 1426.27219764 49 708.796279081 50 1629.74208491 51 434.901189434 52 1367.8300756 53 800.328947676 54 1602.7841183 55 466.650977551 56 1426.27230322 57 708.796118075 58 1629.74208668 59 434.901187396 60 1367.83007161 61 800.32895407 62 1602.78411471 63 466.650981886 64 1426.2723107 65 708.79610667 66 1629.74208681 67 434.901187252 68 1367.83007133 69 800.328954523 70 1602.78411445 71 466.650982193 72 1426.27231123 73 708.796105862 74 1629.74208682 75 434.901187242 76 1367.83007131 77 800.328954555 78 1602.78411443 79 466.650982215 80 1426.27231127 81 708.796105804 82 1629.74208682 83 434.901187241 84 1367.83007131 85 800.328954557 86 1602.78411443 87 466.650982216 88 1426.27231127 89 708.7961058 90 1629.74208682 91 434.901187241 92 1367.83007131 93 800.328954557 94 1602.78411443 95 466.650982216 96 1426.27231127 97 708.7961058 98 1629.74208682 99 434.901187241 100 1367.83007131 |
# The values do not converge upon a specific number. They move around wildly, but there does appear to be some sortof pattern to the values.
|
|
print"Project 1, E"
print "Using algebra to find the equilibrium points, we solved the equation"
print "x=4.2x*e^(-(x/1000)^2)and found that when r=4.2, the two equilibrium are"
print "0 and 1000*ln(4.2)^(5/21)";
Project 1, E Using algebra to find the equilibrium points, we solved the equation x=4.2x*e^(-(x/1000)^2)and found that when r=4.2, the two equilibrium are 0 and 1000*ln(4.2)^(5/21) Project 1, E Using algebra to find the equilibrium points, we solved the equation x=4.2x*e^(-(x/1000)^2)and found that when r=4.2, the two equilibrium are 0 and 1000*ln(4.2)^(5/21) |
p=var('p')
f(p)=diff(4.2*p*e^(-(p/1000)^2), p);
print(abs(f(0))-1,"unstable");print(abs(f(1000*ln(4.2)^(5/21)))-1,"unstable");
(3.20000000000000, 'unstable') (0.761428551041400, 'unstable') (3.20000000000000, 'unstable') (0.761428551041400, 'unstable') |
cobweb(4.2*x*e^(-(x/1000)^2),x,700,10, 1000);
print" x=1000*ln(4.2)^(5/21) is an unstable equilibrium point";
cobweb(4.2*x*e^(-(x/1000)^2),x,.1,10, 1000);
print" x=0 is an unstable equilibrium point";
|
x=1000*ln(4.2)^(5/21) is an unstable equilibrium point
x=0 is an unstable equilibrium pointprint " n p"
P = vector (RDF,10000)
P[0] = 200.
for n in range(1,1001): P[n]=4.2*P[n-1]*e^(-(P[n-1]/1000)^2);
for n in range(0,1001): print "%3s %11s" %(n,P[n])
WARNING: Output truncated! full_output.txt n p 0 200.0 1 807.063128888 2 1767.17149708 3 326.792348175 4 1233.5065361 5 1131.34490506 6 1321.21311039 7 968.536392983 8 1592.09786825 9 530.127075422 10 1681.03768599 11 418.368814788 12 1475.00319612 13 703.357665266 14 1801.25343148 15 294.950826679 16 1135.57815553 17 1313.49127471 18 982.665878794 19 1571.39880271 20 558.643647046 21 1717.30861892 22 377.831836828 23 1375.78100005 24 870.517632045 25 1713.63408723 26 381.806651303 27 1386.06286846 28 852.468951809 29 1731.10982616 30 363.166459356 31 1336.82936405 32 940.138984022 33 1631.49317978 34 478.456580206 35 1598.35346096 36 521.693385181 37 1669.03445532 38 432.425077422 39 1506.4364257 40 654.084127368 41 1790.93741613 42 304.332895937 43 1165.13025712 44 1259.08957301 45 1083.46785208 46 1406.84387638 47 816.460816515 48 1760.67965671 49 333.134680385 50 1252.1942155 51 1096.3555653 52 1384.14214318 53 855.829255425 54 1727.98574366 55 366.449782975 56 1345.68785274 57 924.145294368 ... 941 380.016326384 942 1381.44638305 943 860.554327937 944 1723.49161859 945 371.210277907 946 1358.39093701 947 901.369011308 948 1679.97327076 949 419.602358293 950 1477.82380562 951 698.85774397 952 1801.05806823 953 295.126461351 954 1136.13660636 955 1312.47110779 956 984.536623635 957 1568.60702555 958 562.561119077 959 1721.77204281 960 373.043432837 961 1363.23792641 962 892.730602624 963 1689.8608293 964 408.240457694 965 1451.39524673 966 741.605162062 967 1797.08876868 968 298.711960597 969 1147.49370546 970 1291.653136 971 1022.89782102 972 1508.93510464 973 650.251250598 974 1789.36603188 975 305.781377256 976 1169.64161372 977 1250.72128083 978 1099.11049583 979 1379.2527345 980 864.406854536 981 1719.74081507 982 375.217149536 983 1368.95306396 984 882.583651689 985 1701.02144027 986 395.675638873 987 1421.00563833 988 792.305933635 989 1776.29218728 990 318.032640495 991 1207.24210545 992 1180.56074188 993 1230.41362164 994 1137.14102042 995 1310.63540363 996 987.905162308 997 1563.55072428 998 569.699071704 999 1729.58326299 1000 364.768181165 WARNING: Output truncated! full_output.txt n p 0 200.0 1 807.063128888 2 1767.17149708 3 326.792348175 4 1233.5065361 5 1131.34490506 6 1321.21311039 7 968.536392983 8 1592.09786825 9 530.127075422 10 1681.03768599 11 418.368814788 12 1475.00319612 13 703.357665266 14 1801.25343148 15 294.950826679 16 1135.57815553 17 1313.49127471 18 982.665878794 19 1571.39880271 20 558.643647046 21 1717.30861892 22 377.831836828 23 1375.78100005 24 870.517632045 25 1713.63408723 26 381.806651303 27 1386.06286846 28 852.468951809 29 1731.10982616 30 363.166459356 31 1336.82936405 32 940.138984022 33 1631.49317978 34 478.456580206 35 1598.35346096 36 521.693385181 37 1669.03445532 38 432.425077422 39 1506.4364257 40 654.084127368 41 1790.93741613 42 304.332895937 43 1165.13025712 44 1259.08957301 45 1083.46785208 46 1406.84387638 47 816.460816515 48 1760.67965671 49 333.134680385 50 1252.1942155 51 1096.3555653 52 1384.14214318 53 855.829255425 54 1727.98574366 55 366.449782975 56 1345.68785274 57 924.145294368 ... 941 380.016326384 942 1381.44638305 943 860.554327937 944 1723.49161859 945 371.210277907 946 1358.39093701 947 901.369011308 948 1679.97327076 949 419.602358293 950 1477.82380562 951 698.85774397 952 1801.05806823 953 295.126461351 954 1136.13660636 955 1312.47110779 956 984.536623635 957 1568.60702555 958 562.561119077 959 1721.77204281 960 373.043432837 961 1363.23792641 962 892.730602624 963 1689.8608293 964 408.240457694 965 1451.39524673 966 741.605162062 967 1797.08876868 968 298.711960597 969 1147.49370546 970 1291.653136 971 1022.89782102 972 1508.93510464 973 650.251250598 974 1789.36603188 975 305.781377256 976 1169.64161372 977 1250.72128083 978 1099.11049583 979 1379.2527345 980 864.406854536 981 1719.74081507 982 375.217149536 983 1368.95306396 984 882.583651689 985 1701.02144027 986 395.675638873 987 1421.00563833 988 792.305933635 989 1776.29218728 990 318.032640495 991 1207.24210545 992 1180.56074188 993 1230.41362164 994 1137.14102042 995 1310.63540363 996 987.905162308 997 1563.55072428 998 569.699071704 999 1729.58326299 1000 364.768181165 |
# The values do not converge upon a specific number. They move around wildly, but there does appear to be some sortof pattern to the values.
|
|
# Problem 2 A
# s will define the size of the mammal in centimeters
s = vector (RDF, 34)
s = (122, 125, 145, 87.5, 135, 135, 21.5, 700, 72.5, 120, 320, 41, 505, 135, 38, 75, 77.5, 130, 57.5, 165, 125, 330, 34, 50.5, 52, 135, 65, 67, 150, 210, 140, 125.5, 76, 130)
# w will define the weight of the mammal in kilograms
w = vector (RDF, 34)
w = (66, 50, 46, 16, 323, 136, 0.7515, 2932.5, 9, 113, 400, 2.67, 1600, 57.5, 1.4, 10.5, 6.1, 90, 7.2, 70, 30, 300, 3.099, 6, 1.007, 48.5, 8, 13, 124.5, 93, 200, 48.5, 20, 325)
|
|
# Problem 2 B
Q = [[s[i], w[i]] for i in range(34)]
Q;
pplot= scatter_plot(Q);
pplot.show()
|
|
# The above graph suggests that the size/weight relationship is probably exponential. The least squares regression proves this...
n=34;
x = [s[i] for i in range(n)];
y = [w[i] for i in range(n)];
xx = [s[i]*s[i] for i in range(n)];
xy = [s[i]*w[i] for i in range(n)];
one = [1 for i in range(n)];
A =matrix(RDF, 2,3, [[sum(xx), sum(x), sum(xy)], [sum(x), sum(one), sum(y)]]);
A
[ 1326096.75 4897.5 3347376.05725] [ 4897.5 34.0 7158.2275] [ 1326096.75 4897.5 3347376.05725] [ 4897.5 34.0 7158.2275] |
AE = A.echelon_form();
AE
[ 1.0 4.33680868994e-19 3.73207163737] [ 0.0 1.0 -327.046863059] [ 1.0 4.33680868994e-19 3.73207163737] [ 0.0 1.0 -327.046863059] |
m= AE[0,2];
b= AE[1,2];
(m,b)
(3.73207163737, -327.046863059) (3.73207163737, -327.046863059) |
Q = [[x[i],y[i]] for i in range(n)]
pplot= scatter_plot(Q);
z = var('z')
leastsquaresline= plot(m*z+b, (z,-3,750),rgbcolor=hue(0.8))
show(pplot+leastsquaresline)
|
|
# It is clear that the Least Squares Regression is not a good fit for the data.
# Problem 2 C
# We will try a log-log plot.
S = [ln(s[i]) for i in range(n)];
W = [ln(w[i]) for i in range(n)];
Q = [[S[i],W[i]] for i in range(n)]
pplot= scatter_plot(Q);
pplot.show()
|
|
# Now we will try to fit a Least Squares Regression to the new logarithmic data.
n=34;
x = [S[i] for i in range(n)];
y = [W[i] for i in range(n)];
xx = [S[i]*S[i] for i in range(n)];
xy = [S[i]*W[i] for i in range(n)];
one = [1 for i in range(n)];
A =matrix(RDF, 2,3, [[sum(xx), sum(x), sum(xy)], [sum(x), sum(one), sum(y)]]);
AE = A.echelon_form();
m= AE[0,2];
b= AE[1,2];
(m,b)
(2.53863936438, -8.31024575707) (2.53863936438, -8.31024575707) |
Q = [[x[i],y[i]] for i in range(n)]
pplot= scatter_plot(Q);
z = var('z')
leastsquaresline= plot(m*z+b, (z,2,7.5),rgbcolor=hue(0.8))
show(pplot+leastsquaresline)
|
|
# It is clear that this line is a much better fit to the data.
|
|
# Problem 2 D
# We can use this log-log regression to find an equation that relates the raw data.
m = 2.53863936438
b = -8.31024575707
s = vector (RDF, 34)
s = (122, 125, 145, 87.5, 135, 135, 21.5, 700, 72.5, 120, 320, 41, 505, 135, 38, 75, 77.5, 130, 57.5, 165, 125, 330, 34, 50.5, 52, 135, 65, 67, 150, 210, 140, 125.5, 76, 130)
w = vector (RDF, 34)
w = (66, 50, 46, 16, 323, 136, 0.7515, 2932.5, 9, 113, 400, 2.67, 1600, 57.5, 1.4, 10.5, 6.1, 90, 7.2, 70, 30, 300, 3.099, 6, 1.007, 48.5, 8, 13, 124.5, 93, 200, 48.5, 20, 325)
q = [[s[i], w[i]] for i in range(34)]
q;
pplot= scatter_plot(q);
z = var('z')
NEWline= plot((e^b)*(z^m), (z,-3,750),rgbcolor=hue(0.8))
show(pplot+NEWline)
|
|
# The new equation is w = (e^b)*(s^m) where m = 2.53863936438 and b = -8.31024575707
|
|
# Problem 2 E
# We can approximate m = 2.53863936438 to be (5/2). This tells us that there is, perhaps, a biological/physical relationship that can predict weight based on size. If we consider the general relationship between volume and length (or height), the weight should increase cubically in proportion to the size. The value we got (2.53863936438) is rather close to 3.
|
|
|
|