diffgame.py (8811B)
1 # Editor's note: Simulations for a problem in my high school differential equations class. 2 # I remember very little about how this works, but it produces some flashy GIFs. 3 # Written around September 2022. 4 # 5 # The Game 6 # 7 # 1. Alice chooses three scalars, $\alpha$, $\beta$, and $\gamma$, as well as two points, 8 # $(x_1, y_1)$ and $(x_2, y_2)$, in the first quadrant. 9 # 2. Bob chooses a function $f(x)$ that satisfies $\lim_{x \rightarrow 0} f(x) = b < \infty$ such 10 # that there exists a solution to $(\alpha + \beta x + \gamma x^2)y' + \lambda y = f(x)$ that passes 11 # through one or both of ${(x_1, y_1), (x_2, y_2)}$. Let the set of points that Bob's solution 12 # passes through be $S$. 13 # 3. Alice finds a distinct solution from Bob's that passes through the points in $S$. 14 # 15 # The last player to make a move wins. That is, Bob can win by Alice failing at step 3, while Alice 16 # can win either by Bob failing at step 2 or by herself succeeding at step 3. 17 # 18 # We want to determine who wins with what strategies for which values of $\lambda$ (set before the 19 # game as a parameter) and design a computational environment to let us experience the fun of 20 # differential games ourselves. 21 22 import matplotlib.pyplot as plt 23 import numpy as np 24 from random import * 25 import math 26 seed(10) 27 from numba import jit 28 from sympy import * 29 from zope.interface import * 30 from typing import Tuple, Callable, NewType, Optional 31 from tqdm import tqdm 32 from enum import Enum 33 from matplotlib.animation import FuncAnimation 34 import sys 35 36 # Utility funcs 37 def euler(y_0, deriv, x_range, step=0.001): 38 eps = 1e-3 39 ys = [y_0] 40 xs = [x_range[0]] 41 # print(xs, l, f(xs[-1]), l*ys[-1]) 42 for i in range(int((x_range[1] - x_range[0]) / step)): 43 # if abs(ys[-1]) == 0: continue 44 ys.append(ys[-1] + deriv(xs[-1], ys[-1]) * step) 45 xs.append(xs[-1] + step) 46 return xs, ys 47 48 49 # globals 50 51 l = 1 52 53 # interfaces 54 55 56 Point = NewType("Point", Tuple[float, float]) 57 58 59 class AliceStrat(Interface): 60 def round1() -> Tuple[float, float, float, Point, Point]: 61 """Picks alpha, beta, gamma constants and the two points solutions can/must go through.""" 62 63 def round3(f: Callable[[float], float]) -> Optional[Callable[[float], float]]: 64 """Picks a unique solution for the differential equation.""" 65 66 67 class BobStrat(Interface): 68 def round2( 69 a: float, b: float, c: float, p1: Point, p2: Point 70 ) -> Optional[Tuple[Callable[[float], float], Callable[[float], float]]]: 71 """Chooses a real-valued function f(x) and finds a solution to the differential equation.""" 72 73 74 # strats 75 76 77 @implementer(AliceStrat) 78 class LipschitzAndPray: 79 def round1(): 80 """Hope that y(x) is not Lipschitz continuous at 0.""" 81 return 0, random(), 0, (0, random()), (0, random()) 82 83 def round3(f): 84 """Pray it doesn't get this far - if it does, give up.""" 85 return None 86 87 88 @implementer(AliceStrat) 89 class NaiveRoots: 90 def round1(): 91 root = random() 92 b, a = uniform(0, root), uniform(root, 1) 93 return root**2, -2*root, 1, (b, random()), (a, random()) 94 95 def round3(f): 96 return None 97 98 @implementer(AliceStrat) 99 class Roots: 100 def round1(): 101 root1, root2 = random(), random() 102 while root1 == root2: # reinit if unlucky 103 root1, root2 = random(), random() 104 return root1*root2, -(root1+root2), 1, (root1, random()), (root2, random()) 105 106 def round3(f): 107 return None 108 109 110 @implementer(AliceStrat) 111 class Random: 112 def round1(): 113 return random(), random(), random(), (random(), random()), (random(), random()) 114 115 def round3(f): 116 return None 117 118 @implementer(AliceStrat) 119 class Rude: 120 """Be mean.""" 121 122 def round1(): 123 big = sys.float_info.max 124 return big, big, big, (big, big), (big, big) 125 126 def round3(f): 127 return sys.float_info.min 128 129 130 @implementer(BobStrat) 131 class PrecomputedAnalyticZero: 132 """Assumes a,b,c are nonzero and picks f(x) = 0 to make things simple. 133 Plugs into a precomputed analytical solution to get y(x). 134 """ 135 def round2(a, b, c, p1, p2): 136 try: 137 scary_term = lambda x: math.exp( 138 -(2 * l * math.atan((b + 2 * c * x) / math.sqrt(4 * a * c - b**2))) 139 / math.sqrt(4 * a * c - b**2) 140 ) 141 k = p1[1] / scary_term(p1[0]) 142 return lambda x: 0, lambda x: k * scary_term(x) 143 except: 144 return None 145 146 147 @implementer(BobStrat) 148 class Analytic: 149 """Uses SymPy to solve the differential equation.""" 150 def __init__(self, f): 151 self.f = f 152 153 def round2(self, a, b, c, p1, p2): 154 y = Function("y") 155 sa, sb, sc, sl, sx = symbols("a b c l x") 156 eq = Eq(Derivative(y(sx), sx), (self.f-sl * y(sx) / (sa + sb * sx + sc * sx**2))) 157 eq = eq.subs({sa: a, sb: b, sc: c, sl: l}) 158 sol = dsolve(eq) 159 print(sol) 160 k = p1[1]/sol.subs({Symbol("C1"): 1, sx: p1[0]}).rhs 161 sol = sol.subs({Symbol("C1"): k}) 162 return lambda x: 0, lambda x: sol.subs(sx, x).rhs 163 164 165 @implementer(BobStrat) 166 class Euler: 167 def __init__(self, f): 168 self.f = f 169 170 """Approximates a solution numerically using Euler's method.""" 171 def round2(self,a,b,c,p1,p2): 172 try: 173 deriv = lambda x, y: (self.f(x) - l * y) / (a + b * x + c * x**2) 174 xs, ys = euler(p1[1], deriv, (p1[0], 1)) 175 except: 176 return None 177 178 # print(ys) 179 def near_analytic(x): 180 for i in range(len(xs))[:-1]: 181 if x >= xs[i] and x < xs[i+1]: 182 return ys[i] 183 return np.nan 184 185 return self.f, near_analytic 186 187 # verifying code 188 189 def verify_diffeq(y, deriv, p1, p2): 190 """Verifies a potential solution. 191 - Checks that it satisfies the differential equation 192 - Checks that it passes through at least one of the points 193 """ 194 # verify it satisfies the diffeq 195 h = 0.001 196 for x in np.arange(p1[0], p2[0], 0.01): # FIXME, super naive 197 if abs((y(x + h) - y(x)) / h - deriv(x, y)) > 0.001: 198 # print((y(x+h) - y(x))/h, deriv(x,y)) 199 return None 200 201 epsilon = 0.001 202 return (abs(y(p1[0])-p1[1]) < epsilon, abs(y(p2[0]) - p2[1]) < epsilon) 203 204 class Result(Enum): 205 WIN = (1,) 206 GIVEUP = (2,) 207 NO_POINT = (3,) 208 INVALID_EQ = (4,) 209 ALICE_WIN = (5,) 210 211 def __str__(self): 212 return self.name 213 214 alice = NaiveRoots 215 bob = Euler(lambda x: 0) 216 217 fig = plt.figure() 218 ax = plt.axes() 219 ax.set_title("NaiveRoots (Alice) vs. Euler (Bob)") 220 ax.set_xlim(-.1, 1.1) 221 ax.set_ylim(-.1, 1.1) 222 line, = ax.plot([], [], lw=2) 223 scatter = ax.scatter([], []) 224 res_text = ax.text(0.05, 0.9, '', transform=ax.transAxes) 225 stats_text = ax.text(0.05, 0.05, '', transform=ax.transAxes) 226 results = [] 227 228 def play_game(i): 229 a, b, c, p1, p2 = alice.round1() 230 if p1[0] <= 0 or p1[1] <= 0 or p2[0] <= 0 or p2[1] <= 0: 231 return # TODO make less bad 232 scatter.set_offsets([p1, p2]) 233 234 bob_choice = bob.round2(a, b, c, p1, p2) 235 if bob_choice != None: 236 f, y = bob_choice 237 238 result = Result.WIN 239 240 xs = np.arange(0, 1, 0.01) 241 ys = list(map(y, xs)) 242 deriv = lambda x, y: (f(x) - l * y(x)) / (a + b * x + c * x**2) 243 pts = verify_diffeq(y, deriv, p1, p2) 244 245 # this is so hacky 246 # for i in ys: 247 # if i.is_real is None: 248 # result = Result.GIVEUP 249 # pts = None 250 # line.set_data([], []) 251 # break 252 # else: 253 line.set_data(xs, ys) 254 if pts is not None: 255 # TODO implement forcing Alice to go through the same points as Bob 256 if True in pts: 257 y2 = alice.round3(f) 258 if y2 is not None and verify_diffeq(y2, deriv, p1, p2): 259 result = Result.ALICE_WIN 260 else: 261 result = Result.NO_POINT 262 else: 263 result = Result.INVALID_EQ 264 else: 265 result = Result.GIVEUP 266 267 results.append(result) 268 res_text.set_text(str(result)) 269 res_text.set_c("g" if str(result) == "WIN" else "r") 270 271 wins = round(results.count(Result.WIN)/len(results) * 100, 2) 272 no_points = round(results.count(Result.NO_POINT)/len(results) * 100, 2) 273 giveups = round(results.count(Result.GIVEUP)/len(results) * 100, 2) 274 invalid_eqs = round(results.count(Result.INVALID_EQ)/len(results) * 100, 2) 275 alice_wins = round(results.count(Result.ALICE_WIN)/len(results) * 100, 2) 276 277 lines = [f"WIN {wins}%", f"GIVEUP {giveups}%", 278 f"INVALID_EQ {invalid_eqs}%", f"ALICE_WIN {alice_wins}%"] 279 stats_text.set_text("\n".join(lines)) 280 281 return [line, scatter, res_text, stats_text] 282 283 anim = FuncAnimation(fig, play_game, frames=20, interval=20, blit=True) 284 anim.save("./roots_vs_euler.gif")