```python
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
```


```python
%matplotlib notebook
```


```python
def average_fecundity(U, T, R, P, S):
    return 2 * ((R + P - T - S) * U**2 + (T + S - 2 * P) * U + P)

def equilibrium_coordinators_share(U, T, R, P, S):
    F = average_fecundity(U, T, R, P, S)
    return 2 * ((R - S) * U**2 + S * U) / F
```


```python
U = sym.symbols('U', real=True, nonnegative=True)
T, R, P, S = sym.symbols('T, R, P, S', real=True, positive=True)
x = equilibrium_coordinators_share(U, T, R, P, S)
```


```python
first_derivative = sym.lambdify((U, T, R, P, S), sym.diff(x, U, 1), modules=["numpy"])
second_derivative = sym.lambdify((U, T, R, P, S), sym.diff(x, U, 2), modules=["numpy"])
```


```python
_, ax = plt.subplots(1,1)
Us = np.logspace(-6, 0, 1000)
payoffs = [6, 2, 4, 1]
xs = equilibrium_coordinators_share(Us, *payoffs)
x_primes = first_derivative(Us, *payoffs)
x_prime_primes = second_derivative(Us, *payoffs)
ax.plot([0,1], [0,1], "k--")
ax.plot(xs, Us)
ax.plot(xs, x_primes, label=r"$\frac{\partial U}{\partial x}$")
ax.plot(xs, x_prime_primes, label=r"$\frac{\partial^2 U}{\partial x^2}$")
ax.set_xlabel(r"$x^*$")
ax.set_ylabel(r"$U(x^*)$", rotation="horizontal")
ax.set_title("T={}, R={}, P={}, S={}".format(*payoffs))
ax.legend()

```


    <IPython.core.display.Javascript object>









    <matplotlib.legend.Legend at 0x14e991e80>




```python
first_derivative(0, *payoffs), first_derivative(1, *payoffs)
```




    (0.5, 1.3333333333333335)




```python
# slope of the L-locus at x=0, U=0
1 / sym.diff(x, U).subs({U: 0})
```




    P/S




```python
# slope of the L-locus at x=1, U=1
1 / sym.together(sym.diff(x, U).subs({U: 1}))
```




    R/T




```python
sym.factor(sym.diff(x, U), U)
```




    -(-P*S + U**2*(2*P*R - P*S - R*T) + U*(-2*P*R + 2*P*S))/(P + U**2*(P + R - S - T) + U*(-2*P + S + T))**2




```python
numerator, denominator = sym.fraction(sym.factor(sym.diff(x, U), U))
```


```python
# always positive!
denominator
```




    (P + U**2*(P + R - S - T) + U*(-2*P + S + T))**2




```python
numerator
```




    P*S - U**2*(2*P*R - P*S - R*T) - U*(-2*P*R + 2*P*S)



Numerator polynomial opens up $\iff PS + RT > 2 PR$ 


```python
numerator_poly = sym.poly(numerator, U)
```


```python
sym.factor(numerator_poly.discriminant())
```




    4*P*R*(P*R - S*T)



Discriminant of the numerator polynomial will be negative (implying polynomial has imaginary roots!) $\iff$ $\frac{P}{S} < \frac{T}{R}$


```python
sym.solve(numerator_poly, U)
```




    [(sqrt(P)*sqrt(R)*sqrt(P*R - S*T) - P*(R - S))/(-2*P*R + P*S + R*T),
     -(sqrt(P)*sqrt(R)*sqrt(P*R - S*T) + P*(R - S))/(-2*P*R + P*S + R*T)]



Sufficient conditions for first derivative to be strictly positive are that numerator polynomial opens up and discriminant of the polynomial is negative.


```python

```


```python
second_derivative(0, *payoffs), second_derivative(1, *payoffs)
```




    (1.5, 0.4444444444444445)



$$ (-U*(S + U*(R - S))*(-2*P + S + T + 2*U*(P + R - S - T)) + (S + 2*U*(R - S))*(P + U**2*(P + R - S - T) + U*(-2*P + S + T))) > 0$$


```python
sym.simplify(sym.diff(x, U, 2).subs({U: 0}))
```




    (2*P*(R - S) - 2*S*(-2*P + S + T))/P**2




```python
sym.simplify(sym.diff(x, U, 2).subs({U: 1}))
```




    (-2*P*R - 2*R*T + 2*S*T + 2*T**2)/R**2




```python
sym.together(sym.diff(x, U, 2))
```




    2*(R*(P + U**2*(P + R - S - T) + U*(-2*P + S + T))**2 - S*(P + U**2*(P + R - S - T) + U*(-2*P + S + T))**2 - U*(S + U*(R - S))*(P + U**2*(P + R - S - T) + U*(-2*P + S + T))*(P + R - S - T) + U*(S + U*(R - S))*(-2*P + S + T + 2*U*(P + R - S - T))**2 - (S + 2*U*(R - S))*(P + U**2*(P + R - S - T) + U*(-2*P + S + T))*(-2*P + S + T + 2*U*(P + R - S - T)))/(P + U**2*(P + R - S - T) + U*(-2*P + S + T))**3




```python
numerator, denominator = sym.fraction(sym.simplify(sym.diff(x, U, 2)))
```


```python
# always positive!
denominator
```




    (P + U**2*(P + R - S - T) + U*(-2*P + S + T))**3




```python
numerator_poly = sym.poly(numerator, U)
```


```python
sym.factor(numerator_poly.discriminant())
```




    432*P**2*R**2*(P*R - S*T)**2*(P + R - S - T)**2*(4*P*R - S**2 - 2*S*T - T**2)




```python
sym.solve(numerator_poly, U)
```




    [-(-3*P*R + 3*P*S)/(3*(2*P*R - P*S - R*T)) - (9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)/(3*(27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)) - (27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)/3,
     -(-3*P*R + 3*P*S)/(3*(2*P*R - P*S - R*T)) - (9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)/(3*(-1/2 - sqrt(3)*I/2)*(27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)) - (-1/2 - sqrt(3)*I/2)*(27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)/3,
     -(-3*P*R + 3*P*S)/(3*(2*P*R - P*S - R*T)) - (9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)/(3*(-1/2 + sqrt(3)*I/2)*(27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)) - (-1/2 + sqrt(3)*I/2)*(27*P*S*(-3*P*R + 3*P*S)/(2*(2*P*R - P*S - R*T)**2) + (-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + sqrt(-4*(9*P*S/(2*P*R - P*S - R*T) + (-3*P*R + 3*P*S)**2/(2*P*R - P*S - R*T)**2)**3 + (27*P*S*(-3*P*R + 3*P*S)/(2*P*R - P*S - R*T)**2 + 2*(-3*P*R + 3*P*S)**3/(2*P*R - P*S - R*T)**3 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2))**2)/2 + 27*(P**2*R + P**2*S - P*S**2 - P*S*T)/(2*(2*P**2*R - P**2*S + 2*P*R**2 - 3*P*R*S - 3*P*R*T + P*S**2 + P*S*T - R**2*T + R*S*T + R*T**2)))**(1/3)/3]




```python

```

b647d163a719bce884e97598469e84dfb4e1b923
ipynb
Jupyter Notebook
Untitled2.ipynb
davidrpugh/evolution-of-social-behavior
52e8370613ab4f246759f28e578a631b86d8e712
MIT
2019-01-16T04:19:13.000Z
2020-01-28T11:45:17.000Z
Untitled2.ipynb
davidrpugh/evolution-of-social-behavior
52e8370613ab4f246759f28e578a631b86d8e712
MIT
2018-01-17T11:09:59.000Z
2018-01-17T11:09:59.000Z
Untitled2.ipynb
davidrpugh/evolution-of-social-behavior
52e8370613ab4f246759f28e578a631b86d8e712
MIT
2020-05-12T13:43:08.000Z
2021-04-19T07:14:14.000Z
Qwen/Qwen-72B
1. YES
2. YES