computer_assignment2
March 8, 2023
[5]: import numpy as np
import matplotlib.pyplot as plt
def find_a(lambda_val):
# define function to find a
def f(a):
return a * np.exp(-a) - lambda_val * np.exp(-lambda_val)
# use bisection method to find root of f between 0 and lambda_val
a_left = 0
a_right = lambda_val
while abs(a_right - a_left) > 1e-6:
a_mid = (a_left + a_right) / 2
if f(a_mid) == 0:
return a_mid
elif f(a_mid) * f(a_left) < 0:
a_right = a_mid
else:
a_left = a_mid
return a_mid
# initialize arrays for lambda and pi values
lambda_vals = np.arange(1.05, 3.05, 0.05)
pi_vals = []
# find a and pi values for each lambda
a_vals = []
for lambda_val in lambda_vals:
a = find_a(lambda_val)
pi = round(a/lambda_val, 3)
pi_vals.append(pi)
a_vals.append(a)
# print table of lambda and pi values
print("lambda\t a_vals\t pi")
for i in range(len(lambda_vals)):
1
print("{:.2f}\t{:.3f}\t{:.3f}".format(lambda_vals[i],a_vals[i], pi_vals[i]))
# plot pi vs lambda
plt.plot(lambda_vals, pi_vals)
plt.xlabel("lambda")
plt.ylabel("extinction probability")
plt.show()
lambda a_vals pi
1.05 0.952 0.906
1.10 0.906 0.824
1.15 0.864 0.751
1.20 0.824 0.686
1.25 0.786 0.629
1.30 0.750 0.577
1.35 0.716 0.531
1.40 0.685 0.489
1.45 0.654 0.451
1.50 0.626 0.417
1.55 0.599 0.386
1.60 0.573 0.358
1.65 0.548 0.332
1.70 0.525 0.309
1.75 0.503 0.287
1.80 0.482 0.268
1.85 0.461 0.249
1.90 0.442 0.233
1.95 0.424 0.217
2.00 0.406 0.203
2.05 0.390 0.190
2.10 0.374 0.178
2.15 0.358 0.167
2.20 0.344 0.156
2.25 0.330 0.147
2.30 0.316 0.138
2.35 0.304 0.129
2.40 0.291 0.121
2.45 0.280 0.114
2.50 0.268 0.107
2.55 0.258 0.101
2.60 0.247 0.095
2.65 0.237 0.090
2.70 0.228 0.084
2.75 0.219 0.080
2.80 0.210 0.075
2.85 0.202 0.071
2.90 0.194 0.067
2
2.95 0.186 0.063
3.00 0.179 0.060
[6]: 0.952*np.exp(-0.952)
[6]: 0.36744183528427854
[7]: 1.05*np.exp(-1.05)
[7]: 0.3674346365667131
[ ]: