Computational Approach to Game Theory

Dr. Gerald H. Thomas, Milwaukee School of Engineering

Join the Conversation #WolframTechConf

Abstract:

A course was developed for Junior and Senior undergraduates at MSOE on “A computational approach to decision process theory”. The course was built around Mathematica notebook segments that provided the core mathematics for the basic game theory and differential geometry concepts of the theory using my recently published book Geometry, Language and Strategy, the Dynamics of Decision Processes, Vol. 2 (for a copy, see www.amazon.com/author/strategy).  The students were a mix of software, computer and electrical engineers. No prior experience with Mathematica, game theory or differential geometry was assumed. Based on this approach the students used the provided notebook segments to compute the Nash equilibria and harmonic and transient harmonic behaviors for numerous examples from game theory. In the process students learned the basics of Mathematica and were able to apply Mathematica’s linear programming to get the game theory equilibria and NDSolve to get the harmonic and transient behaviors, learning and applying the concepts of decision process theory; without Mathematica, these would have been formidable programming tasks. This talk will provide details of the approach, samples of the notebook segments and provide lessons learned from teaching the course.

Course Goals [10 week course]

In week 1, I presented the expected outcomes for the course:

Be proficient in Computational Engineering using Mathematica

Be able to understand how the dynamic thinking approach of engineering and physics can be applied to economic problems

Be able to understand and apply the Systems Thinking approach to science and social problems

Be able to understand the basic concepts of the Theory of Games as an equilibrium theory

Go beyond Systems Thinking and understand the extension of the Theory of Games to a dynamic theory

Be able to apply the dynamic theory to Robinson Crusoe dynamics [and single person games]

Be able to apply the dynamic theory to the Prisoner's Dilemma [and two person games]

Be able to understand the more advanced concepts of Game Theory and how they are articulated in a dynamic theory

Be able to look for dynamic behaviors of the theory in data, such as the use of recurrence plots of stock market data

Labs and Homework--How we hoped to achieve the expected outcomes

Lab 1: Basic Hello World Problems: Arithmetic, Lists, Plots

Hm 2: Game Theory Example

Lab 3: Lists and Images: Strings, Text, Sounds, Arrays

Lab 4: Solvers: Solve, RSolve, DSolve

Hm 5: Game Theory Examples continued [using provided notebooks as a base]

Hm 6: none

Hm 7: Decision Process Theory Harmonics [using provided notebooks]

Lab 8-10: (student chosen) Research Topic and Presentation

Basic Approach: lab/homework/final project

For Mathematica, use Wolfram, An Elementary Introduction to the Wolfram Language

For the Game Theory and Decision Process Theory use Thomas, Geometry, Language and Strategy--the Dynamics of Decision Processes, Vol. 2 (for a copy see www.amazon.com/author/strategy) with the following segments as seen in this notebook:

Initializations,

Game Matrix

and Harmonic Flows

Submit all work as a Mathematica Notebook

Midterms and Final were also submitted as email as a Mathematica Notebook from each student

Introduction of Mathematica

Challenges

Fear of calculus

Mixture of backgrounds

Software Engineers were strong in algorithms with a focus on programming languages such as Java, C++ and other higher level languages

Electrical Engineers were strong on mathematical tools that simulate directly the hardware such as Spice and NI MultiSim

Computer Engineers were more focused on embedded systems, and low level languages such as C, Assembler languages along with some use of Spice, NI MultiSim and MatLab

Approach

Use Wolfram’s An Elementary Introduction to the Wolfram Language as a base

Use Mathematica’s Natural Language Support aspects as applied to the current subject matter.

For example TIC-TAC-TOE as a game theory game can be found using Mathematica Wolfram Alpha, giving not only the extensive but intensive description

Goal was to bring all students whichever their background quickly to the point where they could use the initializations and start working on the consequences of the game theory

As students skill progressed, they started to understand the Mathematica language

From the beginning, worked on building the calculus computational skills with the help of NDSolve (students did grasp the concepts of calculus, but were fuzzy on how to use them)

About half way through the course, the students were able to go beyond simple game theory and start to use and understand the notebooks to compute the harmonic behaviors

Students then moved away from decisions as a study of static equilibrium to decisions as a reflection of dynamic processes

At some point, students started to grasp the fundamental constructs of the theory

By the time of the final project, students started stretching their (and my) understanding of how the theory works on practical problems

So by the end of the course I believe that at least some of the students understood and were able to use everything that was in the notebook; all the students could use the key aspects of the notebooks

Deliverables

The midterm and final were done online and submitted by email to me at the end of the exam as Mathematica Notebooks

The labs, homework and final project were all submitted as Mathematica Notebooks.

Part of the grading rubric was that the notebooks used clear formatting and organization. Some of the students became very good at that.

Partial List of Final Projects Created by Students

“Chicken” and other variants: Hawk Dove, ...

Tragedy of the Commons

Variations on Work Wealth Wisdom Model

Baseball strategy Model

Football Strategy Model

Internet Card Trading Game strategies

Dating Game

Initializations

This notebook will use the following formulas. The section is put at the beginning to make sure it is available when needed.

Finite Element Initialization

Needs["NDSolve`FEM`"]

Useful Methods

gToSymG[x_,y_,g_]:=ArrayFlatten[{{0,-Transpose[g],List/@y},{g,0,-List/@x},{-List@y,List@x,0}}]
hedge1Fn[g_]:=Table[1,{k,Length[g]}]
hedge2Fn[g_]:=Table[1,{k,Length[Transpose[g]]}]
hedgeSym[g_]:=Join[hedge1Fn[g],hedge2Fn[g],{1}]
toPos[g_]:=g-Min[g]+1
solP2[g_]:=LinearProgramming[hedge2Fn[g],toPos[-g],hedge1Fn[g]]
solP1[g_]:=LinearProgramming[hedge1Fn[g],toPos@Transpose@g,hedge2Fn[g]]
normStrategy[s_]:=If[Tr[s]≠0,(s/(Tr[s])),"NA"]
gameValueP2[g_]:=((-1)/(solP2[g].hedge2Fn[g]))-Min[-g]+1
gameValueP1[g_]:=(1/(solP1[g].hedge1Fn[g]))+Min[g]-1
symMatrix[g_]:=gToSymG[hedge1Fn[g],hedge2Fn[g],g]
symSol[g_]:=LinearProgramming[hedgeSym[g],toPos[-gToSymG[hedge1Fn[g],hedge2Fn[g],toPos@g]],hedgeSym[g]]
symStrategy1[g_]:=normStrategy@symSol[g][[-Length@hedge1Fn@g-1;;-2]]
symStrategy2[g_]:=normStrategy@symSol[g][[1;;Length@hedge2Fn@g]]
symStrategy1[g_,sol_]:=normStrategy@sol[[-Length@hedge1Fn@g-1;;-2]]
symStrategy2[g_,sol_]:=normStrategy@sol[[1;;Length@hedge2Fn@g]]
pValuation2[g_]:=Transpose[g].normStrategy@solP1[g]
pValuation1[g_]:=g.normStrategy@solP2[g]
dptGame[g_,eScale_]:=gToSymG[pValuation1[g]/eScale,pValuation2[g]/eScale,g]
dptGameWWW[g_,f0_,f1_,f2_]:=gToSymG[f2 pValuation1[g]/f0,f1 pValuation2[g]/f0,g]
dropTimePayoff[f_]:=#[[;;-2]]&/@f[[;;-2]]
numberColumns[g_]:=Length@First@g
numberRows[g_]:=numberColumns@Transpose@g
playerOne[g_]:=DiagonalMatrix@Join[Table[0,{k,numberColumns[g]}],Table[1,{k,numberRows[g]}],{0}]
playerTwo[g_]:=DiagonalMatrix@Join[Table[1,{k,numberColumns[g]}],Table[0,{k,numberRows[g]}],{0}]
gab[g_]:=DiagonalMatrix@Join[Table[-1,{k,numberRows[g]+numberColumns[g]}],{1}]
timeStrategy[g_]:=Diagonal[1-playerOne@g-playerTwo@g]
kthPlayerOneStrategy[k_,g_]:=Join[Table[0,{j,numberColumns@g}],Table[If[j==k,1,0],{j,numberRows@g}],{0}]
kthPlayerTwoStrategy[k_,g_]:=Join[Table[If[j==k,1,0],{j,numberColumns@g}],Table[0,{j,numberRows@g}],{0}]

Setup for differential equations and plots

numActiveStrategies:=numberColumns@gameMatrix+numberRows@gameMatrix
numStrategies:=numActiveStrategies+2
strategyString:=Join[strategyString2,strategyString1]
strategyColors:=Join[colors2,colors1]
allStrategyString:=Join[strategyString,{"time"}]
allStrategyColors:=Join[strategyColors,{Purple}]
myDC[s1_,s2_]:=s1==0||s2==0
hi1[s1_,s2_]:=e1+0e1 Exp[-(((((s1+s2)/(2size))))^2)]
hi2[s1_,s2_]:=e2 + 0Sin[((2Pi s1)/size)]
ha1[s1_,s2_]:=0Sin[((2Pi s2)/size)]
ha2[s1_,s2_]:=0Sin[((2Pi s1)/size)]
ha[s1_,s2_]:=Diagonal[one ha1[s1,s2]+two ha2[s1,s2]]
hap[s1_,s2_]:=Diagonal[one D[ha1[s1,s2],s1]+two D[ha2[s1,s2],s2]]
strategyString1:=Table[StringJoin["P1S",ToString[k]],{k,numberRows[gameMatrix]}]
colors1:=Table[RGBColor[0,((k-1)/(numberRows[gameMatrix]-1)),1-((k-1)/(numberRows[gameMatrix]-1))],{k,numberRows[gameMatrix]}]

strategyString2:=Table[StringJoin["P2S",ToString[k]],{k,numberColumns[gameMatrix]}]
colors2:=Table[RGBColor[1,((k-1)/(numberColumns[gameMatrix])),0],{k,numberColumns[gameMatrix]}]

setups for the differential equations for the harmonic solutions

(*yInit[s1_,s2_]:=initdisruptiveInfluenceFlow(Sqrt[1+(e1^2)+(e2^2)])((1/2)s1+(1/2)s2)+ha[s1,s2]*)
yInit[s1_,s2_]:=(initdisruptiveInfluenceFlow.one s1+initdisruptiveInfluenceFlow.two s2+(metric+one+two).initdisruptiveInfluenceFlow((1/2)s1+(1/2)s2))(Sqrt[1+(e1^2)+(e2^2)])+ha[s1,s2]
initCovariantFlow:=metric.initdisruptiveInfluenceFlow(Sqrt[1+(e1^2)+(e2^2)])+hap[s1,s2]
(* the alternate choice is to assume that we have the harmonic flow along the boundary based on the One Owner Model *)
(*initCovariantFlow:=(1/2)metric.one.flowSolution[s1/2]+(1/2)metric.two.flowSolution[s2/2]+hap[s1,s2]
initCovariantFlow:=metric.flowSolution[(1/2)s1+(1/2)s2]+hap[s1,s2]*)
flow[s_]:=Table[(Subscript[x, k])'[s],{k,numActiveStrategies+1}]
var:=Table[(Subscript[x, k])[s],{k,numActiveStrategies+1}]
covariantV[s_]:=metric.flow[s]
flowEquations:=Table[covariantV'[s][[k]]==((e1 dptGameMatrix1+e2 dptGameMatrix2).flow[s])[[k]],{k,numActiveStrategies+1}]
boundaryEquations:=Join[Table[flow[0][[k]]==initdisruptiveInfluenceFlow[[k]](Sqrt[1+(e1^2)+(e2^2)]),{k,numActiveStrategies+1}],Table[(Subscript[x, k])[0]==0,{k,numActiveStrategies+1}]]

Pressure Model

(*effortModel[x_]:=2-Exp[- (x^2)]*)
effortModel[x_]:=Exp[- (x^2)]
ϕ1[pref_]:=setPressureOne effortModel[setInversePressureRangeOne Tr[(  ((fTwo one)/(Tr[one]))- ((fOne two)/(Tr[two]))).pref]]
ϕ2[pref_]:=setPressureTwo effortModel[setInversePressureRangeTwo Tr[( ((fOne two)/(Tr[two]))- ((fTwo one)/(Tr[one]))).pref]]
gradϕ1[pref_]:=Table[D[ϕ1[pref],(Subscript[y, k])],{k,numActiveStrategies+1}]
gradϕ2[pref_]:=Table[D[ϕ2[pref],(Subscript[y, k])],{k,numActiveStrategies+1}]

Frame Model

frameModel[x_]:=-setLockedBehaviorPotentialWell (1/(1+(((setLockedBehaviorRange x))^2)))
ν[pref_]:=frameModel[Tr[(  ((fTwo one)/(Tr[one]))- ((fOne two)/(Tr[two]))).pref]]
Φ[pref_]:=-(1/(numStrategies-2))ν[pref]
gradν[pref_]:=Table[D[ν[pref],(Subscript[y, k])],{k,numActiveStrategies+1}]
gradΦ[pref_]:=Table[D[Φ[pref],(Subscript[y, k])],{k,numActiveStrategies+1}]

Seasonal Game Model: (Subscript[f, υυ']) and (Subscript[(G^k), υυ'])

characteristicPayoff[g_]:=gToSymG[pValuation1[g]0,pValuation2[g]0,g]
seasonalVectorPotential[pref_]:=-(1/2)characteristicPayoff[characteristicGameMatrix].pref
dSeasonalVectorPotential[pref_]:=Table[Table[D[seasonalVectorPotential[pref][[k]],(Subscript[y, j])]-D[seasonalVectorPotential[pref][[j]],(Subscript[y, k])],{k,numActiveStrategies+1}],{j,numActiveStrategies+1}]
seasonalPayoff[pref_]:=dSeasonalVectorPotential[pref]

Harmonic Perturbations

harmonic[pref_,ω_,var_]:=Laplacian[var@@pref,pref]+(ω^2) var@@pref==0
nDVariables:=Rule@@@Transpose[{prefyS,ToExpression@Table[StringJoin["y",ToString[k]],{k,Length[prefyS]}]}]
region:=Table[{(Subscript[y, k]),0,ySize},{k,Length[nDVariables]}]/.nDVariables

Sample Games

We list here a number of games that you might find of interest. You can feel free to augment this list as needed. Several of the games are taken from (Williams, 1966).

baseGame:={{7,1,3,0,2},{0,1,6,4,2},{1,2,0,5,5}}
rspGame:={{0,1,-1},{-1,0,1},{1,-1,0}}

rspwGame:={{0,1,-1,1/2},{-1,0,1,1/2},{1,-1,0,1/2},{1,1,1,0}}

rspwFairGame:={{0,1,-1,-1/2},{-1,0,1,-1/2},{1,-1,0,-1/2},{1/2,1/2,1/2,0}}

morra:={{0,2,2,-3,0,0,-4,0,0},{-2,0,0,0,3,3,-4,0,0},{-2,0,0,-3,0,0,0,4,4},{3,0,3,0,-4,0,0,-5,0},{0,-3,0,4,0,4,0,-5,0},{0,-3,0,0,-4,0,5,0,5},{4,4,0,0,0,-5,0,0,-6},{0,0,-4,5,5,0,0,0,-6},{0,0,-4,0,0,-5,6,6,0}}

theDate:={{(1/3),0,-(1/3)},{0,(2/3),(1/3)},{(1/3),-(1/3),1}}

colonelBlotto:={{4,2,1,0},{1,3,0,-1},{-2,2,2,-2},{-1,0,3,1},{0,1,2,4}}

catMouse:={{0,0,0,1,0,0,0,1},{0,1,1,1,1,1,1,0},{0,1,1,1,1,1,1,0},{1,1,1,1,1,1,1,0},{0,1,1,1,1,1,1,1},{0,1,1,1,1,1,1,0},{0,1,1,1,1,1,1,0},{1,0,0,0,1,0,0,0}}

p208:={{5,1,3},{2,4,1}}

p102:={{2900,3000,3200},{18000,6000,-2000},{2000,7000,12000}}

p80:={{-4,-2,3,4},{6,5,0,1}}

prisonersDilemma:={{-(1/10),-1},{0,-(9/10)}}
contracts:=prisonersDilemma+1;
PDMafia[x_]:=10prisonersDilemma+{{0,x},{-x,0}}
hornsGame :={{0,2,2,0,0,0,0,0,0},{0,0,0,0,3,3,0,0,0},{0,0,0,0,0,0,0,4,4},{3,0,3,0,0,0,0,0,0},{0,0,0,4,0,4,0,0,0},{0,0,0,0,0,0,5,0,5},{4,4,0,0,0,0,0,0,0},{0,0,0,5,5,0,0,0,0},{0,0,0,0,0,0,6,6,0}}

attackDefenseBlue:={{4,1},{3,4}}
attackDefenseRed:=attackDefenseBlue
hm5Ex1:={{4,1},{3,4}}
hm5Ex2:={{0,-1},{-3,0}}
hm5Ex3:={{250,-150},{-150,350}}
hm5Ex4:={{0,1,-1},{-1,0,1},{1,-1,0}}
hm5Ex5:=morra
wwwModel:={{0,-(1/10)},{-(9/10),0}}
drinkingGame:= {{55,10},{10,110}}
chicken:={{0,-1},{1,-10}} (* row labels are "swerve" and "go straight" *)
chickenSwitch:={{0,11},{9,10}} (*chicken with row labels "go straight" and "swerve" corresponding to econChicken labels *)
econChicken1:= {{0,4},{1,3}} (* row labels are "crash" and "swerve" *)
econChicken2:=-Transpose[econChicken1]
politicalChicken:={{-4,1},{2,0}}
politicalChicken2:={{-1,2},{5,1}}
politicalChicken3:={{-1,2},{10,1}}
dareChicken:={{0,7},{2,6}}  (* row labels are "dare" and "chicken out" *)
dareChickenOut:={{0,4},{1,3}}  (* row labels are "dare" and "chicken out" *)
dareChickenOutExtreme:={{0,4},{1,10}}  (* row labels are "dare" and "chicken out" *)
myChicken:={{-10,4},{3,1}}  (* row labels are "crash" and "swerve" *)
toyChicken:={{-10,-3},{-8,-4}}   (* row labels are "crash" and "swerve" *)
stagHare:={{3,0},{2,1}} (* row labels are "stag" and "hare" *)
stagHarev2:={{3,0},{1,2}} (* row labels are "stag" and "hare" *)
tragedyOfCommons:={{100, 70}, {140, 80}} (*fishing example, 1 fish, 2 fish options *)
DotaSoloPro:={{72,100,100},{96,50,100},{96,100,50}} (*model by Ziwei Chen CS498: stay in line, go left rune, go right rune *)
DotaSoloNoob:={{24,100,100},{96,50,100},{96,100,50}} (*model by Ziwei Chen CS498 *)
tragedyOfCommons:={{100, 70}, {140, 80}} (*Ryan Forward: fishing with one boat, two boat *)
modernMeta := (({{0, -.05, .07, -.06, -.05, -.08}, {.05, 0, -.15, .14, .06, -.09}, {-.07, .15, 0, -.12, -.17, .19}, {.06, -.14, .12, 0, .09, -.11}, {.05, -.06, .17, -.09, 0, -.02}, {.08, .09, -.19, .11, .02, 0}})) (*Brett Peavler Trading Card Game *)

Sample Characteristic Games

unitValue[g_]:=g 0+1

Ownership Model Mechanics Notebook [background material]

Preliminaries on Aspects of Calculus and Notation

From our discussions, it is clear we weren’t familiar with some of the mathematics. Here is a quick refresher course in the relevant calculus.

There is the chain rule from calculus: ((df(x(s)))/dt)=(df/dx)(dx/ds)

The same rule works with two variables: ((df(x(s),y(s)))/ds)=((∂f)/(∂x))(dx/ds)+((∂f)/(∂y))(dy/ds)

Now assume that x and y are coordinates and their derivative along s is the flow vector V: ((df(x(s),y(s)))/ds)=((∂f)/(∂x))(V^x)+((∂f)/(∂y))(V^y)

Now write this in compact notation: ((df(x(s),y(s)))/ds)=(V^a)((∂f)/(∂(x^a)))=(V^a)(∂_a)f

This is the notation used in this notebook.

From today, we observed that the harmonics we are computing correspond to resonant harmonics of the system. In addition there can be forced harmonics, which are not treated in this course (see the textbook however if you are interested)

Abstract

New principle for decision making: Decision Flow moves along the shortest path

The geodesic

Essential Details from Differential Geometry

We argue that decision processes have aspects that are not just deterministic, but share a property of other physical systems in that the behavior may be chaotic in the technical sense: solutions to the equations may be sensitive to initial conditions and produce novel structures such as turbulence. Ordinarily we think of determinism as being related to cause and effect; something that happens now impacts behaviors in the future. An equally important concept is that of elasticity: things that happen nearby will impact the behavior “here” in the future. A simple example is hold a rubber band between your hands. When you move one hand, the other will move, though not at exactly the same instant. Neighbor events impact each other and are thus connected! 

We believe this is also true in decision making. The behavior of neighboring stocks (in the same sector say) will influence each other. This is an extension of cause and effect, which we term elasticity. This aspect is captured in decision process theory, which has all of these characteristics.

 In this notebook we use numerical models to highlight some of these aspects. Since the equations are non-linear and not particularly transparent, it is useful to find simple models that produce the results we expect from the full theory. One such way of looking at the theory is to characterize it as producing behavior that follows equations similar to the motion of charged particles in an electromagnetic field. This way of looking at the theory is based on simplifications that are discussed in Vol. 2.

We found that the total acceleration is determined by the ownership forces that are constraints outside of the geometry:

		μ(Subscript[(Overscript[V, ˙]), μ])=((Subscript[h, μ])^σ)(Subscript[p, σν;λ])(g^(λν))

The ownership forces are on the right. The left hand side is geometric and can be expanded in a particular reference frame to display market and cooperative forces:

		μ(Subscript[(Overscript[V, ˙]), a])=μ(V^b)(∂_b)(Subscript[V, a])−(1/2)μ(∂_a)(Subscript[g, bc])(V^b)(V^c)−μ(Subscript[V, k])(Subscript[F, ab]^k)(V^b)−(1/2)μ(∂_a)(Subscript[γ, kj])(V^k)(V^j)

Frame dependent rotational effects: gradient forces (centripetal and gravitational forces) and velocity dependent forces (Coriolis forces) come from −(1/2)μ(∂_a)(Subscript[g, bc])(V^b)(V^c). To the extent that the metric valuations are constants along the streamlines, we ignore these “gravitational” effects.

The ownership forces on the right are documented elsewhere (see the website:

http://decisionprocesstheory.com/white-papers/vol-2-resources/stationary-ownership-model/

We can consider metric contributions that are consistent with the covariant gauge:

					({{(∂_(υ'))(((Overscript[g, ¯])^(υ'υ))(Sqrt[(Overscript[g, ¯])γ(Subscript[(Overscript[g, ¯]), οο])]))=0}, {n=(Subscript[n, a])+(Subscript[n, i])}, {−2Φ+nΦ+ν=0⇒Φ=−(ν/(n−2))}, {(∂_ο)ν=0}})

We work in the time orthogonal frame, which leads us to the conclusion that the payoffs will in general depend on the preference positions

					({{(Subscript[((Overscript[F, ¯])^k), υυ'])=(Subscript[(F^k), υυ'])+(Subscript[a, υ])(Subscript[(F^k), υ'ο])−(Subscript[a, υ'])(Subscript[(F^k), υο])≡(Subscript[(F^k), υυ'])+(Subscript[(G^k), υυ'])}, {(Subscript[((Overscript[F, ¯])^k), ου])=(Subscript[(F^k), ου])}, {(Subscript[f, υυ'])=(∂_υ)(Subscript[a, υ'])−(∂_(υ'))(Subscript[a, υ])}})

We take our formulas from this site, Sec. 1.8. The results are expressed in terms of the vector fields defined for each player ((Subscript[S, J])^υ)≡(g^(υυ[J]))(Subscript[V, υ[J]]) . So for example, the inactive flows are:

			(μ+(Subscript[φ, J]))((∂(Subscript[V, j[J]]))/(∂(Subscript[s, J])))+(∑_(K≠J))μ((∂(Subscript[V, j[J]]))/(∂(Subscript[s, K])))=−(Subscript[V, j[J]])((∂(Subscript[φ, J]))/(∂(Subscript[s, J])))+(Subscript[V, j[J]])(Subscript[θ, νλ])(p^(νλ))

The parameters (Subscript[s, J]) are the coordinate measures defined along the streamlines  for the player vector field ((Subscript[S, J])^υ). The active flows are similarly determined:

	(∑_K)μ((∂(Subscript[V, υ[J]]))/(∂(Subscript[s, K])))+(Subscript[φ, J])((∂(Subscript[(Overscript[V, ¯]), υ[J]]))/(∂(Subscript[s, J])))=(({{(∂_(υ[J]))(Subscript[φ, J])−(Subscript[(Overscript[V, ¯]), υ[J]])((Overscript[V, ¯])^(υ'[J]))(∂_(υ'[J]))(Subscript[φ, J])+(Subscript[(Overscript[V, ¯]), υ[J]])(Subscript[(Overscript[θ, ¯]), νλ])((Overscript[p, ¯])^(νλ))}, {+(μ+(Subscript[φ, J]))(1−((Overscript[V, ¯])^ο)(Subscript[(Overscript[V, ¯]), ο]))(∂_(υ[J]))Φ}, {+(μ(Subscript[(Overscript[V, ¯]), ο])((Overscript[V, ¯])^ο)+(Subscript[φ, J])(1−(Overscript[ϕ, ¯])(Subscript[(Overscript[ψ, ¯]), J])(1−(Subscript[(Overscript[ψ, ¯]), J]))))(∂_(υ[J]))ν}, {+(∑_K)(Subscript[(Overscript[V, ¯]), ο])(Overscript[ϕ, ¯])((μ+(Subscript[φ, K]))(1−(Subscript[(Overscript[ψ, ¯]), K]))−(∑_(L≠K))μ(Subscript[(Overscript[ψ, ¯]), L]))(Subscript[(Overscript[f, ¯]), υ[J]υ'[K]])((Overscript[V, ¯])^(υ'[K]))}, {+(∑_K)((μ+(Subscript[φ, K]))(Subscript[(Overscript[V, ¯]), k[K]])((Overscript[F, ¯])^(k[K]))(Subscript[ , υ[J]υ'[K]])+(∑_(L≠K))μ(Subscript[V, l[L]])((Overscript[F, ¯])^(l[L]))(Subscript[ , υ[J]υ'[K]]))((Overscript[V, ¯])^(υ'[K]))}, {−(∑_K)(Overscript[ϕ, ¯])((Overscript[V, ¯])^ο)((μ+(Subscript[φ, K]))(1−(Subscript[(Overscript[ψ, ¯]), K]))−(∑_(L≠K))μ(Subscript[(Overscript[ψ, ¯]), L]))(Subscript[V, k[K]])((Overscript[F, ¯])^(k[K]))(Subscript[ , ου[J]])}})).

These equations reflect only the transverse conservation laws of the theory. The final equation is determined from the longitudinal conservation law:

		(∑_K)(((∂μ)/(∂(Subscript[s, K])))+2μ((∂Φ)/(∂(Subscript[s, K]))))=(∑_K)(({{−(Overscript[ϕ, ¯])((Overscript[V, ¯])^ο)(Subscript[φ, K])(Subscript[(Overscript[V, ¯]), k[K]])(∑_(L≠K))((Subscript[(Overscript[ψ, ¯]), L])((Overscript[F, ¯])^(k[K]))(Subscript[ , ου[K]])((Overscript[V, ¯])^(υ[K]))+(1−(Subscript[(Overscript[ψ, ¯]), K]))((Overscript[F, ¯])^(k[K]))(Subscript[ , ου[L]])((Overscript[V, ¯])^(υ[L])))}, {−(Subscript[φ, K])((Overscript[V, ¯])^(υ'[K]))(Subscript[(Overscript[V, ¯]), k[K]])((Overscript[F, ¯])^(k[K]))(Subscript[ , υ'[K]υ])((Overscript[V, ¯])^υ)−(Overscript[ϕ, ¯])(Subscript[(Overscript[V, ¯]), ο])(Subscript[φ, K])(1−(Subscript[(Overscript[ψ, ¯]), K]))((Overscript[V, ¯])^(υ'[K]))(Subscript[f, υ'[K]υ])((Overscript[V, ¯])^υ)}, {+(Subscript[φ, K])((Overscript[V, ¯])^(k[K]))((∂(Subscript[(Overscript[V, ¯]), k[K]]))/(∂(Subscript[s, K])))+(Subscript[φ, K])((Overscript[V, ¯])^(υ[K]))((∂(Subscript[(Overscript[V, ¯]), υ[K]]))/(∂(Subscript[s, K])))−(Overscript[ϕ, ¯])(Subscript[φ, K])(1−(Subscript[(Overscript[ψ, ¯]), K]))(((Overscript[V, ¯])^υ)((∂(Subscript[(Overscript[V, ¯]), υ]))/(∂(Subscript[s, K])))+((Overscript[V, ¯])^k)((∂(Subscript[(Overscript[V, ¯]), k]))/(∂(Subscript[s, K]))))}, {+(Subscript[φ, K])(Subscript[(Overscript[ψ, ¯]), K])((∂ν)/(∂(Subscript[s, K])))−(Overscript[ϕ, ¯])(Subscript[φ, K])(Subscript[(Overscript[ψ, ¯]), K])(1−(Subscript[(Overscript[ψ, ¯]), K]))((Overscript[V, ¯])^υ)(∂_υ)ν−(Subscript[φ, K])(n[K]−(Subscript[(Overscript[ψ, ¯]), K]))((Overscript[V, ¯])^υ)(∂_υ)Φ}, {+(Subscript[φ, K])(n−2(Subscript[(Overscript[ψ, ¯]), K]))((∂Φ)/(∂(Subscript[s, K])))+(Overscript[ϕ, ¯])(Subscript[φ, K])(1−(Subscript[(Overscript[ψ, ¯]), K]))(1−((Overscript[V, ¯])^ο)(Subscript[(Overscript[V, ¯]), ο]))((∂Φ)/(∂(Subscript[s, K])))}})).
		
These expressions form a complete set of equations if we take as given the various player stresses, (Subscript[φ, K]). 

In this notebook we will proceed as follows. We will first define the game matrix to be used for the numerical model; then we will isolate the terms in the above equations that are analogous to electromagnetic theory and look at their numerical solution, the harmonic flows; we will then set up a recursive way to solve the coupled equations above and then we will provide the numerical solution.

Game Matrix--Two Players

In this notebook, we treat any of the games in the initialization section subject to a number of simplifications. We assume that along the shortest path, the metric is constant, that the stresses are constant and that there are only two players. Further details can be found on the website referenced above. As necessary, we comment on the more detailed assumptions in the relevant text below.

Commentary on the Decision Process

Before specifying the game matrix for each player, in this section, the student should describe the decision process. These are the requirements, the objectives for the process under consideration. For example, for this talk the game is economic chicken, where each player has two options: Crash or Swerve [It is optional to redefine the strategyString1 and strategyString2 to meet your specific model; default names are provided].

strategyString1={"OneAimsToCrash","OneSwerves"}
strategyString2={"TwoAimsToCrash","TwoSwerves"}

The assumed payoff matrix for each player is

MatrixForm[econChicken1]

The idea is to start the model off with a disruptive influence, set below.

Choose Game Matrix for each Player

We need to pick a game matrix for each player. The two game matrices will agree on the number of strategies for each player, but may differ otherwise.  The various functions referred to here are defined in the initialization section at the beginning of the notebook. We also define below the projection operator that defines which strategies each player owns: see one and two below.

zeroSumGame=False
MatrixForm[gameMatrix= econChicken1]
MatrixForm[gameMatrix1=gameMatrix]
MatrixForm[gameMatrix2=If[zeroSumGame,gameMatrix,-Transpose[gameMatrix]]] (*for non-zero sum games--another option is that there may be a different gameMatrix *)
MatrixForm[metric=gab[gameMatrix]]
MatrixForm[one=playerOne[gameMatrix]]
MatrixForm[two=playerTwo[gameMatrix]]
valueScale=4 (* sets scale of player valuations *)
sMax=10(* max "invariant" path length that includes time: (ds^2)=(dt^2)- (∑_k)(dy^k)(dy^k) *)
fOne=(10/10) (*population weight of player 1*)
fTwo=((10 )/10)(*population weight of player 2 *)
μ0=5(* energy density needed in ownership model *)
e1=(1/2) (* player 1 engagement, contributes to composite payoff, use the sign corresponding to the gradient capacitance model *)
e2=If[zeroSumGame,e1,-e1]  (* player 2 engagement, contributes to composite payoff, we use the gradient capacitance model *)
MatrixForm[characteristicGameMatrix=(1/sMax) unitValue[gameMatrix]]
setLockedBehaviorRange=3(4/sMax)
setLockedBehaviorPotentialWell=(0/10)
setPressureOne=(1/2)(1/(Tr[one]))
setPressureTwo=(1/2)(1/(Tr[two]))
setInversePressureRangeOne=0(4/sMax)
setInversePressureRangeTwo=0(4/sMax)
setMaxPreferenceForHarmonicCalculation=sMax
setForcedHarmonicTimePeriod=(14/4)sMax
setForcedHarmonicFrequency=((2Pi)/setForcedHarmonicTimePeriod)
meshFactor=5 (* mesh=meshFactor((size^2)/1000) is used for the FEM calculation for the ownership model*)

Direct Method Evaluation of gameMatrix strategies and values (from player 1’s frame of reference)

ps1=normStrategy@solP1@gameMatrix1
ps2=normStrategy@solP2@gameMatrix1
gameValueP1@gameMatrix1//N
gameMatrix1.ps2
ps1.gameMatrix1

ps1=normStrategy@solP1@gameMatrix2
ps2=normStrategy@solP2@gameMatrix2
gameValueP1@gameMatrix2//N
gameMatrix2.ps2
ps1.gameMatrix1

Player One--Game Matrix

Decision process theory uses the symmetrized form of the game, which is accomplished here for player one. If there is more than one null vector, you may need to be careful about which one you pick.

MatrixForm[dptGameMatrix1=dptGameWWW[gameMatrix1,valueScale,fOne,fTwo]]
nullSpace1=NullSpace@dptGameMatrix1

Player Two--Game Matrix

Here we have the symmetrized form for player two. Again be careful about which null vector you pick if there is more than one.

MatrixForm[dptGameMatrix2=dptGameWWW[gameMatrix2,valueScale,fOne,fTwo]]
nullSpace2=First@NullSpace@dptGameMatrix2

Disruptive Influence Flow

We are interested not just in the flow at equilibrium but the flow that starts away from equilibrium; the disruptive influence flow

offset=1
disruptiveInfluenceFlow=offset((kthPlayerOneStrategy[1,gameMatrix])/valueScale)+.5offset((kthPlayerOneStrategy[2,gameMatrix])/valueScale)+0offset((kthPlayerTwoStrategy[1,gameMatrix])/valueScale)+offset((kthPlayerTwoStrategy[2,gameMatrix])/valueScale)+5timeStrategy@gameMatrix 
disruptiveInfluenceFlow.metric.disruptiveInfluenceFlow
(*disruptiveInfluenceFlow=First@nullSpace1*)
initdisruptiveInfluenceFlow=(disruptiveInfluenceFlow/(√(disruptiveInfluenceFlow.metric.disruptiveInfluenceFlow)))//N

Composite Null  Vector

The composite null vector need not be time like; if it is space like the solutions can grow exponentially! However the disturbed flow vector will generate only timelike solutions since the invariant length is conserved.

MatrixForm[e1 dptGameMatrix1+e2 dptGameMatrix2]
Eigenvalues@(metric.(e1 dptGameMatrix1+e2 dptGameMatrix2))//N
NullSpace@(e1 dptGameMatrix1+e2 dptGameMatrix2)
compositeNullVector=First@NullSpace@(e1 dptGameMatrix1+e2 dptGameMatrix2)
compositeNullVector.metric.compositeNullVector

Harmonic Flows--One Owner

Basic Harmonic Results follow from Composite Payoff

If one player owns all of the strategies, then there is only a single payoff matrix to consider (take it to be player one) and the essential equations are:

		(V^b)(∂_b)(Subscript[V, a])=(Subscript[e, 1])(Subscript[(F^1), ab])(V^b)+(Subscript[e, 2])(Subscript[(F^2), ab])(V^b).

We take this form and set the constants to unity. We then have to use the relationship between vectors with upper indices and lower indices, and express the result in terms of a single streamline 

		({{(V^a)=(g^(ab))(Subscript[V, b])}, {(V^a)(∂_a)=(∂/(∂s))}}).

This provides the streamline equation for harmonic flow in decision process theory:

		(∂/(∂s))((Subscript[g, ab])(V^b))=(Subscript[e, 1])(Subscript[(F^1), ab])(V^b)+(Subscript[e, 2])(Subscript[(F^2), ab])(V^b).

The nice thing about mathematics is that it expresses relationships without regard to the domain of interest, relationships that may in fact be similar or identical to some other region of interest. In this case the equation above is the equation for a charged particle in a magnetic and electric field, expressed however in what may be an unfamiliar way.  In any case we can solve this quite easily, which we do below. We set the engagement for each player here as well.

xSolution[s_]=NDSolveValue[Join[flowEquations,boundaryEquations],var,{s,0,sMax},MaxSteps->Infinity]
flowSolution[s_]=xSolution'[s];
contraFlowSolution[s_]=metric.flowSolution[s];
Plot[Evaluate[xSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],PlotLabel->"preferences"]
Plot[Evaluate[flowSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],PlotLabel->"flows"]

Feel free to change the initial conditions. If you pick the null vector, you are likely to get just straight lines. You see that just like a charged particle in a magnetic field, these solutions oscillate around the null vector; in electromagnetic theory, the null vector is the direction of the magnetic field. Here are the pairs of flows for each player:

Strategic Efforts (log plot)

The total effort x expended on a strategy is given by the exponential of the preference. Thus a negative preference corresponds to low effort and a positive preference a higher preference. This can be seen in the LogPlof.

LogPlot[Evaluate[Table[Exp[xSolution[s][[k]]],{k,numActiveStrategies}]],{s,0,sMax},AxesLabel->{"s","(y^a)"},PlotLegends->strategyString,PlotStyle->strategyColors,PlotLabel->"strategic efforts"]

Strategic Efforts (linear plot)

Plot[Evaluate[Table[Exp[xSolution[s][[k]]],{k,numActiveStrategies}]],{s,0,sMax},AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],PlotLabel->"strategic efforts"]

Beta behaviors

Initially, the focus was on null vectors that would simultaneously be the flow vectors. In that case, the null vectors needed to be timelike. In general however, we can take any value for the valueScale and the composite null vector can be timelike, spacelike or indeed may not exist. There will be eigenvalues of the payoff matrix however in all cases. Moreover, in all cases we can compute a physically meaningful timelike flow. We might expect however that when there is a null vector, it will correspond to a fixed point, which indicate a stable or unstable equilibrium. We ascribe the word harmonic to a stable equilibrium and another term, possibly transient flow, to an unstable fixed point. Since the flow vector is timelike, we will have a velocity “ beta”, which is the ratio of the flows divided by the time flow, for each strategic direction which is less than unity, and a preference behavior divide by time preference behavior which is bounded. This is despite the fact that both the flows and preferences may grow exponentially.

Plot[Evaluate[flowSolution[s]/flowSolution[s][[-1]]],{s,0,sMax},AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],PlotLabel->"beta"]

Plot[Evaluate[xSolution[s]/xSolution[s][[-1]]],{s,0,sMax},AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],PlotLabel->"preferences/time"]

Helix Behaviors in phase space plots

For example we can, with an appropriate choice of initial conditions, see evidence of a helix-like behavior in the phase space plots:

TableForm[Table[ParametricPlot[{xSolution[s][[k]],flowSolution[s][[j]]},{s,0,sMax},PlotStyle->Table[allStrategyColors[[k]],{i,numActiveStrategies-k}],PlotLabel->"(y^a) vs (V^b)",AxesLabel->{StringJoin["y",strategyString[[k]]],StringJoin["V",strategyString[[j]]]},AspectRatio->1],{k,numActiveStrategies},{j,numActiveStrategies}]]

We may see the same effect in the scatter plot of one flow against another:

TableForm[Table[ParametricPlot[{flowSolution[s][[k]],flowSolution[s][[j]]},{s,0,sMax},PlotStyle->Table[allStrategyColors[[k]],{i,numActiveStrategies-k}],PlotLabel->"(V^a) vs (V^b)",AxesLabel->{strategyString[[k]],strategyString[[j]]},AspectRatio->1],{k,numActiveStrategies-1},{j,k+1,numActiveStrategies}]]

Preferences

We can do the same thing with the preferences

TableForm[Table[ParametricPlot[{xSolution[s][[k]],xSolution[s][[j]]},{s,0,sMax},PlotStyle->Table[allStrategyColors[[k]],{i,numActiveStrategies-k}],PlotLabel->"(y^a) vs (y^b)",ImageSize->Medium,AxesLabel->{strategyString[[k]],strategyString[[j]]},AspectRatio->1

],{k,numActiveStrategies-1},{j,k+1,numActiveStrategies}]]

Player efforts: relative and total

Tr[one]

Plot[Tr[((one/(Tr[one]))-(two/(Tr[two]))).flowSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotStyle->{Thick,Red},PlotLabel->"relative player effort flow"]
Plot[Tr[(((one+two)/(Tr[one+two]))).flowSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotStyle->{Thick,Blue},PlotLabel->"total player effort flow"]

Plot[Tr[((one/(Tr[one]))-(two/(Tr[two]))).xSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotStyle->{Thick,Red},PlotLabel->"relative player effort"]
Plot[Tr[(((one+two)/(Tr[one+two]))).xSolution[s]],{s,0,sMax},AxesLabel->Automatic,PlotStyle->{Thick,Blue},PlotLabel->"total player effort"]

Player Strategy Differences

Flows

TableForm[Table[Plot[flowSolution[s][[k]]-flowSolution[s][[j]],{s,0,sMax},PlotStyle->Table[allStrategyColors[[k]],{i,numActiveStrategies-k}],AxesLabel->{s,StringJoin["V",strategyString[[k]],"-","V",strategyString[[j]]]},AspectRatio->1],{k,numActiveStrategies-1},{j,k+1,numActiveStrategies}]]

Preferences

TableForm[Table[Plot[xSolution[s][[k]]-xSolution[s][[j]],{s,0,sMax},PlotStyle->Table[allStrategyColors[[k]],{i,numActiveStrategies-k}],AxesLabel->{s,StringJoin["y",strategyString[[k]],"-","y",strategyString[[j]]]},AspectRatio->1],{k,numActiveStrategies-1},{j,k+1,numActiveStrategies}]]

Interpretation of results

[Put in this section the interpretation of the results on harmonic flow. The sections abvove might help to highlight your conclusions.]

Advanced Material

Not all material was covered in depth

Hints as to further directions were given also in notebook segments

For example the ownership properties of the energy momentum tensor were introduced and briefly discussed

Ownership Model and Setup for Linear Recursion

Ownership

For two players, we now have two independent streamlines that we can follow, corresponding to moving along the strategies owned by player one or by player two. The equations form a non-linear partial differential equations, which currently Mathematica can’t deal with directly. What we can do is make an approximation that linearizes the equations by setting the non-linear terms to initial guesses for the unknown functions. We then solve the linearized partial differential equation and obtain the next guess. It may not be perfect and we can’t guarantee it will always work. When it does work we then have a solution.

This section can be run once and makes the initial guess. The next section can then be run multiple times. It is set up to reset the value of the guess at the end. We have found convergent solutions. This section however is a work in progress, with the hope of finding an effective and quick way to get a full solution. We know that Mathematica developers are also working on a more permanent and less ad hoc method.

We set the active flows to the player one null flow values. We can experiment with these values below. We take nominal values for the player engagement and again suggests playing with these parameters. The pressure is no longer a constant so we have to compute

		μ((Subscript[s, 1]),(Subscript[s, 2]))

We take an initial value for this to be large to start with, since that minimizes the effect of changes to the source and may help obtain convergent solutions. This may be a problem currently if there are any algebraic errors in the formulas. As these get worked out, other choices may work just as effectively.

(*lear[s1,s2]*)
sub=Join[Table[(Subscript[u, k])[s1,s2]->initCovariantFlow[[k]],{k,numActiveStrategies}],Table[(Subscript[yOwner, k])[s1,s2]-> yInit[s1,s2][[k]],{k,numActiveStrategies}],{uI1[s1,s2]->hi1[s1,s2],uI2[s1,s2]->hi2[s1,s2],μ[s1,s2]->μ0}]
θpHatOld[s1_,s2_]=0;
numberIterations=0
v0Track[s1_,s2_]=1

Mesh and Region

The functions depend on the distances ((Subscript[s, 1]),(Subscript[s, 2])) along the respective streamlines and reflect how each player makes his or her play. The shape of the region of interest may play a role here, just as such shapes play a role in electromagnetic statics problems. We have started with a square and will note when we come to the results, that the shape influences the behavior of the solutions as much if not at times more than the equations. We suggest playing with different shapes, such as the commented out parallelogram.

On a technical note, the partial differential equations can be solved using the Finite Element Method, which explains the initialization step at the beginning of the chapter and the reference and setup of that method below. One can play with the mesh size using MaxCellMeasure. If you choose a mesh too small, you will run out of memory and may hang your computer. If you make the mesh too big the resultant curves will be erratic and spikey, indicating very poor convergence. Mathematica makes some reasonable choices if you choose not to use this method. It may be mandatory however if you go away from a simple square or rectangle.

It is worth noting that the value of the fields on the boundary are set here by myDC, standing for my Dirichlet Condition. This is a Boolean variable that you set to be true on that portion of the boundary you think provides the initial condition. The equations appear not to be Elliptical and so can be specified by a choice on only a portion of that boundary. This is a fun item to change and play with.

size=sMax
mesh=meshFactor((size^2)/1000)
(*Ω1=Parallelogram[{0,0},{{size,size},{-size,size}}]*)
Ω1=Rectangle[{0,0},{size,size}]
Ω=ToElementMesh[Ω1,MaxCellMeasure->mesh]
Ω["Wireframe"]

(*yInit[s1_,s2_]=one.xSolution[s1/2]+two.xSolution[s2/2]+ha[s1,s2]+(metric+one+two).(xSolution[s1/2]+xSolution[s2/2])*)
(*yInit[s1_,s2_]=xSolution[s1/2]+xSolution[s2/2]+ha[s1,s2]*)
(*yInit[s1_,s2_]=xSolution[(1/2)s1+(1/2)s2]+ha[s1,s2]*)
boundaryCondition=Join[Table[(Subscript[u, k])[s1,s2]==initCovariantFlow[[k]] ,{k,numActiveStrategies}],{uI1[s1,s2]==hi1[s1,s2],uI2[s1,s2]== hi2[s1,s2],μ[s1,s2]==μ0 },Table[(Subscript[yOwner, k])[s1,s2]==yInit[s1,s2][[k]],{k,numActiveStrategies}],{(Subscript[yOwner, numActiveStrategies+1])[s1,s2]==yInit[s1,s2][[numActiveStrategies+1]]}]

Harmonic Flows--Two Owners

Iterations

To solve non-linear PDE, we linearize the PDE by starting with trial values of the unknown factors, then iterate to get a better approximation of these values to obtain our solutions, starting with setting the mesh and boundary. For those doing iterations, we are now on the following iteration:

numberIterations

Pressure Model

prefy=Table[(Subscript[y, k]),{k,numActiveStrategies+1}]
ϕ1[prefy]
ϕ2[prefy]
gradϕ1[prefy]
gradϕ2[prefy]

Here are the pressure and gradients

preferenceList=Table[(Subscript[yOwner, k])[s1,s2],{k,numActiveStrategies+1}]

preferenceListP=preferenceList/.sub;
ϕ1P[s1_,s2_]=Simplify[ϕ1[preferenceListP]];
ϕ2P[s1_,s2_]=Simplify[ϕ2[preferenceListP]];
gradϕ1P[s1_,s2_]=gradϕ1[prefy]/.Table[(Subscript[y, k])->preferenceListP[[k]],{k,numActiveStrategies+1}];
gradϕ2P[s1_,s2_]=gradϕ2[prefy]/.Table[(Subscript[y, k])->preferenceListP[[k]],{k,numActiveStrategies+1}];
gradϕP[s1_,s2_]=one.gradϕ1P[s1,s2]+two.gradϕ2P[s1,s2];

grads1ϕ1[s1_,s2_]=gradϕ1P[s1,s2].flowOneParameter[s1,s2];
grads2ϕ2[s1_,s2_]=gradϕ2P[s1,s2].flowTwoParameter[s1,s2];

For each strategy direction, there is a pressure and coordinate curve gradient associated with that

pressureP[s1_,s2_]=one.Table[1,{k,numActiveStrategies+1}]ϕ1P[s1,s2]+two.Table[1,{k,numActiveStrategies+1}]ϕ2P[s1,s2];

pressurePGrad[s1_,s2_]=one.Table[1,{k,numActiveStrategies+1}]grads1ϕ1[s1,s2]+two.Table[1,{k,numActiveStrategies+1}]grads2ϕ2[s1,s2];

Similarly, for each strategy

u[s1_,s2_]:=Tr[preferenceListP.one-preferenceListP.two]

Frame Model

ν[prefy]
Φ[prefy]
gradν[prefy]
gradΦ[prefy]

νP[s1_,s2_]=Simplify[ν[preferenceListP]];
ΦP[s1_,s2_]=Simplify[Φ[preferenceListP]];
gradνP[s1_,s2_]=gradν[prefy]/.Table[(Subscript[y, k])->preferenceListP[[k]],{k,numActiveStrategies+1}];
gradΦP[s1_,s2_]=gradΦ[prefy]/.Table[(Subscript[y, k])->preferenceListP[[k]],{k,numActiveStrategies+1}];

Variables

We have to list out all the variables that appear as well as any helper functions that may be needed later. You will have one flow variable for every active strategy, called ownership flow, and one for each of the two players representing their engagement flow. You also have the variable for the energy density μ. The variables and derivations are given in the web source cited in the introduction. We will also need to know the values of these variables as determined at the last iteration, which are labeled as Parameter, Parameters or simply with a “P”. Variables in Mathematica that will be used in the PDE will not have values but, like in algebra, are unknowns that are used in the equations.

uI1Parameter[s1_,s2_]=uI1[s1,s2]/.sub;
uI2Parameter[s1_,s2_]=uI2[s1,s2]/.sub;
μP[s1_,s2_]=μ[s1,s2]/.sub;

ownershipVariables:=Table[(Subscript[u, k])[s1,s2],{k,numActiveStrategies}]
inactiveOwners:={uI1[s1,s2],uI2[s1,s2]}
owner[s1_,s2_]=Join[ownershipVariables,{0}];
ownerParameters[s1_,s2_]=owner[s1,s2]/.sub;
ψ1[s1_,s2_]:=(-((uI1Parameter[s1,s2])^2)-ownerParameters[s1,s2].one.ownerParameters[s1,s2])Exp[-2ΦP[s1,s2]]
ψ2[s1_,s2_]:=(-((uI2Parameter[s1,s2])^2)-ownerParameters[s1,s2].two.ownerParameters[s1,s2])Exp[-2ΦP[s1,s2]]
v0[s1_,s2_]:=Exp[νP[s1,s2]]Sqrt[1-ψ1[s1,s2]-ψ2[s1,s2]] (*covariant form *)
v0Old[s1_,s2_]=v0[s1,s2];
ownershipFlow[s1_,s2_]:=Join[ownershipVariables,{v0[s1,s2]}]
flowOne[s1_,s2_]:=Exp[-2ΦP[s1,s2]]Inverse[metric].one.ownershipFlow[s1,s2] (*contravariant form for flow one and two *)
flowTwo[s1_,s2_]:=Exp[-2ΦP[s1,s2]]Inverse[metric].two.ownershipFlow[s1,s2]
flowOneParameter[s1_,s2_]=flowOne[s1,s2]/.sub;
flowTwoParameter[s1_,s2_]=flowTwo[s1,s2]/.sub;

sJ[s1_,s2_]=Diagonal[playerOne[gameMatrix]s1+playerTwo[gameMatrix]s2]
sK[s1_,s2_]=Diagonal[playerOne[gameMatrix]s2+playerTwo[gameMatrix]s1]

Add preferences to the list

variableList=Join[ownershipVariables,preferenceList,{μ[s1,s2]},inactiveOwners]

Test this definition for where it will be used

Table[D[(Subscript[u, k])[s1,s2],sJ[s1,s2][[k]]],{k,numActiveStrategies}]
Table[D[(Subscript[u, k])[s1,s2],sK[s1,s2][[k]]],{k,numActiveStrategies}]

D[(Subscript[s, 1]) (Subscript[s, 1]),(Subscript[s, 1])]

Not sure why I needed to do this, but this forces v0 to be evaluated with the current version of “sub”

Plot3D[v0Old[s1,s2],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},PlotRange->Full,ImageSize->Large,PlotLabel->"Covariant Flow (Subscript[V, ο])Old((Subscript[s, 1]),(Subscript[s, 2]))"]

Pressure Plots

Plot3D[ϕ1P[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","ϕ1"},PlotRange->Full]
Plot3D[ϕ2P[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","ϕ2"},PlotRange->Full]
Plot3D[μP[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","μ"},PlotRange->Full]

Frame Plots

Plot3D[νP[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","ν"},PlotRange->Full]
Plot3D[ΦP[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","Φ"},PlotRange->Full]
Plot3D[μP[s1,s2],{s1,0,size},{s2,0,size},ImageSize->Large,AxesLabel->{"s1","s2","μ"},PlotRange->Full]

Seasonal Payoff Field

MatrixForm[dSeasonalVectorPotential[prefy]]
MatrixForm[seasonalPayoff[prefy]]
Simplify[seasonalVectorPotential[prefy]]
Simplify[seasonalVectorPotential[prefy].one]
Simplify[seasonalVectorPotential[prefy].two]

seasonalVectorPotentialP[s1_,s2_]=seasonalVectorPotential[preferenceListP];
gaugeCorrection1[s1_,s2_]=Table[seasonalVectorPotentialP[s1,s2][[k]]dptGameMatrix1[[j]][[-1]]-seasonalVectorPotentialP[s1,s2][[j]]dptGameMatrix1[[k]][[-1]],{k,numActiveStrategies+1},{j,numActiveStrategies+1}];
gaugeCorrection2[s1_,s2_]=Table[seasonalVectorPotentialP[s1,s2][[k]]dptGameMatrix2[[j]][[-1]]-seasonalVectorPotentialP[s1,s2][[j]]dptGameMatrix2[[k]][[-1]],{k,numActiveStrategies+1},{j,numActiveStrategies+1}];

seasonalVectorPotentialP[(size/5),(size/5)]
MatrixForm[gaugeCorrection1[(size/5),4(size/5)]]
MatrixForm[gaugeCorrection2[(size/5),(size/5)]]

backGroundω=Max@Abs@Eigenvalues@characteristicGameMatrix

seasonalPayoffP[s1_,s2_]=seasonalPayoff[prefy]/.Table[(Subscript[y, k])->preferenceListP[[k]],{k,numActiveStrategies+1}];
MatrixForm[seasonalPayoffP[0,0]]
bEigen=Eigenvalues[seasonalPayoffP[0,0]]
If[backGroundω==0,∞,Chop[((2.Pi I)/(bEigen[[1]]))]]

Plot3D[u[s1,s2],{s1,0,sMax},{s2,0,sMax},AxesLabel->Automatic,PlotLabel->"u(s1,s2)"]
Plot3D[gaugeCorrection1[s1,s2][[3]][[1]],{s1,0,sMax},{s2,0,sMax},AxesLabel->Automatic,PlotLabel->"((G^1))_31(s1,s2)"]
Plot3D[gaugeCorrection2[s1,s2][[3]][[1]],{s1,0,sMax},{s2,0,sMax},AxesLabel->Automatic,PlotLabel->"((G^2))_31(s1,s2)"]

MatrixForm[efff1P[s1_,s2_]=((v0[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))((μP[s1,s2]+ϕ1P[s1,s2])(1-ψ1[s1,s2]) - μP[s1,s2] ψ2[s1,s2] ) seasonalPayoffP[s1,s2]];
MatrixForm[efff2P[s1_,s2_]=((v0[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))((μP[s1,s2]+ϕ2P[s1,s2])(1-ψ2[s1,s2])- μP[s1,s2] ψ1[s1,s2]  )seasonalPayoffP[s1,s2]];
SeasonalPayoffContribution[s1_,s2_]:=efff1P[s1,s2].flowOne[s1,s2]+efff2P[s1,s2].flowTwo[s1,s2]

Functions that are “known” when doing the PDE including Payoffs

As we go through the iterations, we will need to set certain functions so that the PDE is a linear equation. These are some of the functions

MatrixForm[effF1P[s1_,s2_]=(μP[s1,s2]+ϕ1P[s1,s2])uI1Parameter[s1,s2] (dptGameMatrix1 +gaugeCorrection1[s1,s2])+ μP[s1,s2] uI2Parameter[s1,s2] (dptGameMatrix2 +gaugeCorrection2[s1,s2])];
MatrixForm[effF2P[s1_,s2_]=(μP[s1,s2]+ϕ2P[s1,s2])uI2Parameter[s1,s2] (dptGameMatrix2 +gaugeCorrection2[s1,s2]) + μP[s1,s2] uI1Parameter[s1,s2]  (dptGameMatrix1 +gaugeCorrection1[s1,s2])];
payoffContribution[s1_,s2_]:=effF1P[s1,s2].flowOne[s1,s2]+effF2P[s1,s2].flowTwo[s1,s2]

subD1=Rule@@@Transpose[{D[variableList,s1],D[variableList/.sub,s1]}][[All,;;2]];
subD2=Rule@@@Transpose[{D[variableList,s2],D[variableList/.sub,s2]}][[All,;;2]];
totalSub=Join[sub,subD1,subD2];

Valuation Field

For the value fields, which are the time components of the player payoff matrices, we will need to write their contribution to the flow equations in terms of the PDE variables.

valueFieldsP=+Table[(((μP[s1,s2]+ϕ1P[s1,s2])(1-ψ1[s1,s2])-μP[s1,s2] ψ2[s1,s2])/(v0[s1,s2]))uI1Parameter[s1,s2] (dptGameMatrix1 +gaugeCorrection1[s1,s2])[[k]][[-1]]+(((μP[s1,s2]+ϕ2P[s1,s2])(1-ψ2[s1,s2])-μP[s1,s2]  ψ1[s1,s2])/(v0[s1,s2]))uI2Parameter[s1,s2](dptGameMatrix2 +gaugeCorrection1[s1,s2])[[k]][[-1]],{k,numActiveStrategies}];
valueFields=valueFieldsP;

Frame Contributions to Active Flow Equations and to Longitudinal Equation

ψJP[s1_,s2_]=ψ1[s1,s2]Table[1,{k,numActiveStrategies+1}].one+ψ2[s1,s2]Table[1,{k,numActiveStrategies+1}].two;

frame=Table[(μP[s1,s2]+pressureP[s1,s2][[k]])(ψ1[s1,s2]+ψ1[s1,s2])gradΦP[s1,s2][[k]] +(μP[s1,s2](1-ψ1[s1,s2]-ψ1[s1,s2])+pressureP[s1,s2][[k]](1-((ψJP[s1,s2][[k]](1-ψJP[s1,s2][[k]]))/(1-ψ1[s1,s2]-ψ2[s1,s2]))))gradνP[s1,s2][[k]],{k,numActiveStrategies}];

rightHandSideν=ϕ1P[s1,s2]ψ1[s1,s2]flowOneParameter[s1,s2].gradνP[s1,s2]+ϕ2P[s1,s2]ψ2[s1,s2]flowTwoParameter[s1,s2].gradνP[s1,s2]-(1/(1-ψ1[s1,s2]-ψ2[s1,s2]))(flowOneParameter[s1,s2]+flowTwoParameter[s1,s2]).gradνP[s1,s2](ϕ1P[s1,s2]ψ1[s1,s2](1-ψ1[s1,s2])+ϕ2P[s1,s2]ψ2[s1,s2](1-ψ2[s1,s2]));
rightHandSideΦ=-(ϕ1P[s1,s2](Tr[one]+1-ψ1[s1,s2])+ϕ2P[s1,s2](Tr[two]+1-ψ2[s1,s2]))(flowOneParameter[s1,s2]+flowTwoParameter[s1,s2]).gradΦP[s1,s2]+ϕ1P[s1,s2](numActiveStrategies+2-2ψ1[s1,s2])flowOneParameter[s1,s2].gradΦP[s1,s2]+ϕ2P[s1,s2](numActiveStrategies+2-2ψ2[s1,s2])flowTwoParameter[s1,s2].gradΦP[s1,s2]+((ψ1[s1,s2]+ψ2[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))(ϕ1P[s1,s2](1-ψ1[s1,s2])flowOneParameter[s1,s2].gradΦP[s1,s2]+ϕ2P[s1,s2](1-ψ2[s1,s2])flowTwoParameter[s1,s2].gradΦP[s1,s2]);
rightHandSidef=-((v0[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))(ϕ1P[s1,s2](1-ψ1[s1,s2])flowOneParameter[s1,s2].seasonalPayoffP[s1,s2]+ϕ2P[s1,s2](1-ψ2[s1,s2])flowTwoParameter[s1,s2].seasonalPayoffP[s1,s2]).(flowOneParameter[s1,s2]+flowTwoParameter[s1,s2]);
long6=-rightHandSideν-rightHandSideΦ;

Longitudinal Field Equations

the longitudinal equations are particularly lengthy and are summarized here. We have experiment with various choices in order to get solutions that converge in a reasonable time. One approach is to take just the value of the contribution based on the last iteration. That currently is the winner in terms of convergence. Another approach might be to expand that solution to the next order in the difference between the unknown values and the last iteration. This might lead to instabilities however.

MatrixForm[(one+two).Inverse[metric]]

matrixDiffP=Table[ϕ1P[s1,s2]uI1Parameter[s1,s2]dptGameMatrix1[[i]][[j]] -ϕ2P[s1,s2]uI2Parameter[s1,s2]dptGameMatrix2[[i]][[j]],{i,numActiveStrategies+1},{j,numActiveStrategies+1}];

The following might be useful to do a test on the contributions to the energy density equation.

longP=Table[0,{k,6}];
longP[[1]]=Exp[-2ΦP[s1,s2]](ϕ1P[s1,s2]uI1Parameter[s1,s2]D[uI1Parameter[s1,s2],s1]+ϕ2P[s1,s2]uI2Parameter[s1,s2]D[uI2Parameter[s1,s2],s2]);
longP[[2]]=Exp[-2ΦP[s1,s2]](ϕ1P[s1,s2]ownerParameters[s1,s2].one.D[ownerParameters[s1,s2],s1]+ϕ2P[s1,s2]ownerParameters[s1,s2].two.D[ownerParameters[s1,s2],s2]);
longP[[3]]=Exp[-2ΦP[s1,s2]](((1-ψ1[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ1P[s1,s2]ownerParameters[s1,s2].(one+two).Inverse[metric].D[ownerParameters[s1,s2],s1]+((1-ψ2[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ2P[s1,s2]ownerParameters[s1,s2].(one+two).Inverse[metric].D[ownerParameters[s1,s2],s2]-((1-ψ1[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ1P[s1,s2](uI1Parameter[s1,s2]D[uI1Parameter[s1,s2],s1]+uI2Parameter[s1,s2]D[uI2Parameter[s1,s2],s1])-((1-ψ2[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ2P[s1,s2](uI1Parameter[s1,s2]D[uI1Parameter[s1,s2],s2]+uI2Parameter[s1,s2]D[uI2Parameter[s1,s2],s2]));
longP[[4]]=flowOneParameter[s1,s2].matrixDiffP.flowTwoParameter[s1,s2];
longP[[5]]=(1/(v0[s1,s2]))(ϕ1P[s1,s2]uI1Parameter[s1,s2]dptGameMatrix1 [[-1]].(ψ2[s1,s2]flowOneParameter[s1,s2]+(1-ψ1[s1,s2])flowTwoParameter[s1,s2])+ϕ2P[s1,s2]uI2Parameter[s1,s2]dptGameMatrix2 [[-1]].(ψ1[s1,s2]flowTwoParameter[s1,s2]+(1-ψ2[s1,s2])flowOneParameter[s1,s2]));

longP[[6]]=-2μP[s1,s2](flowOneParameter[s1,s2]+flowTwoParameter[s1,s2]).gradΦP[s1,s2]+long6;

We do verification checks here, particularly helpful if we put in equilibrium values for the active flows

test[s1_,s2_]=longP[[1]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[1]]",PlotRange->Full]
test[s1_,s2_]=longP[[2]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[2]]",PlotRange->Full]
test[s1_,s2_]=longP[[3]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[3]]",PlotRange->Full]
test[s1_,s2_]=longP[[2]]+longP[[3]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[2]]+longP[[3]]",PlotRange->Full]
test[s1_,s2_]=longP[[4]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[4]]",PlotRange->Full]
test[s1_,s2_]=longP[[5]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[5]]",PlotRange->Full]
test[s1_,s2_]=longP[[6]];
Plot3D[test[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"longP[[6]]",PlotRange->Full]

This is the form that we use however assuming metric is Minkowski:

long=Table[0,{k,6}];
long[[1]]=Exp[-2ΦP[s1,s2]](ϕ1P[s1,s2]uI1Parameter[s1,s2]D[uI1[s1,s2],s1]+ϕ2P[s1,s2]uI2Parameter[s1,s2]D[uI2[s1,s2],s2]);
long[[2]]=Exp[-2ΦP[s1,s2]](ϕ1P[s1,s2]ownerParameters[s1,s2].one.D[owner[s1,s2],s1]+ϕ2P[s1,s2]ownerParameters[s1,s2].two.D[owner[s1,s2],s2]);
long[[3]]=Exp[-2ΦP[s1,s2]](((1-ψ1[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ1P[s1,s2]ownerParameters[s1,s2].(one+two).Inverse[metric].D[owner[s1,s2],s1]+((1-ψ2[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ2P[s1,s2]ownerParameters[s1,s2].(one+two).Inverse[metric].D[owner[s1,s2],s2]-((1-ψ1[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ1P[s1,s2](uI1Parameter[s1,s2]D[uI1[s1,s2],s1]+uI2Parameter[s1,s2]D[uI2[s1,s2],s1])-((1-ψ2[s1,s2])/(1-ψ1[s1,s2]-ψ2[s1,s2]))ϕ2P[s1,s2](uI1Parameter[s1,s2]D[uI1[s1,s2],s2]+uI2Parameter[s1,s2]D[uI2[s1,s2],s2]));
long[[4]]=(1/2)flowOneParameter[s1,s2].matrixDiffP.flowTwo[s1,s2]+(1/2)flowOne[s1,s2].matrixDiffP.flowTwoParameter[s1,s2];
long[[5]]=(1/(v0[s1,s2]))(ϕ1P[s1,s2]uI1[s1,s2]dptGameMatrix1 [[-1]].(ψ2[s1,s2]flowOneParameter[s1,s2]+(1-ψ1[s1,s2])flowTwoParameter[s1,s2])+ϕ2P[s1,s2]uI2[s1,s2]dptGameMatrix2 [[-1]].(ψ1[s1,s2]flowTwoParameter[s1,s2]+(1-ψ2[s1,s2])flowOneParameter[s1,s2]));

long[[6]]=-2μ[s1,s2](flowOneParameter[s1,s2]+flowTwoParameter[s1,s2]).gradΦP[s1,s2]+long6;

We have set up the helper functions so that the source for the energy density can be expressed in terms of the appropriate fields. We have literally taken the equation and written each line as a contribution and then summed them up below.

θpHatP[s1_,s2_]=Sum[longP[[k]],{k,6}];
Plot3D[θpHatP[s1,s2],{s1,0,sMax},{s2,0,sMax},PlotLabel->"θpHatP]",PlotRange->Full]
θpHat[s1_,s2_]=Sum[long[[k]],{k,6}];

Energy Density Field Equations

The energy density equation (reduced form) is:

μPDE={D[μ[s1,s2],s1]+D[μ[s1,s2],s2]+θpHat[s1,s2]==0};

Inactive Field Equations

We also have all the equations needed to write the inactive flow equations:

inactivePDE={(μP[s1,s2]+ϕ1P[s1,s2])D[uI1[s1,s2],s1]+μP[s1,s2]  D[uI1[s1,s2],s2]-θpHat[s1,s2]  uI1Parameter[s1,s2]+grads1ϕ1[s1,s2]uI1[s1,s2]==0,(μP[s1,s2]+ϕ2P[s1,s2])D[uI2[s1,s2],s2]+ μP[s1,s2]  D[uI2[s1,s2],s1]-θpHat[s1,s2] uI2Parameter[s1,s2]+grads2ϕ2[s1,s2]uI2[s1,s2]==0};

For the single ownership model, the inactive flows were constant along the streamline. This is no longer a given.

Preference Equations

prefCondition=(Table[D[(Subscript[yOwner, k])[s1,s2],s1]==flowOne[s1,s2][[k]],{k,numActiveStrategies+1}].one+Table[D[(Subscript[yOwner, k])[s1,s2],s2]==flowTwo[s1,s2][[k]],{k,numActiveStrategies+1}].two)[[1;;numActiveStrategies]];
timePDE={D[(Subscript[yOwner, numActiveStrategies+1])[s1,s2],s1]+D[(Subscript[yOwner, numActiveStrategies+1])[s1,s2],s2]==v0[s1,s2]};

preferencePDE=Join[prefCondition[[1;;numActiveStrategies]],timePDE];

Streamline Solutions

We have done all of the above preliminaries so that the actual PDE for Mathematica can be written in a fairly condensed form. You may want to be careful about changes you make elsewhere since they may cause this equation to fail. By experience, the most likely cause is that you have introduced a syntax error and something is no longer defined. Look at the error messages for some indication of what the solver thinks is missing.

pdeEqns=Join[Table[(μP[s1,s2]+pressureP[s1,s2][[k]])D[(Subscript[u, k])[s1,s2],sJ[s1,s2][[k]]]+μP[s1,s2]  D[(Subscript[u, k])[s1,s2],sK[s1,s2][[k]]]-θpHat[s1,s2]  ownerParameters[s1,s2][[k]]+pressurePGrad[s1,s2][[k]](Subscript[u, k])[s1,s2]==gradϕP[s1,s2][[k]]+payoffContribution[s1,s2][[k]]+SeasonalPayoffContribution[s1,s2][[k]]+valueFields[[k]]+frame[[k]],{k,numActiveStrategies}],preferencePDE,inactivePDE,μPDE];

streamlines[s1_,s2_]=NDSolveValue[{pdeEqns,DirichletCondition[boundaryCondition,myDC[s1,s2]]},variableList,{s1,s2}∈Ω,MaxSteps->∞]

Strategic Flows, Energy Density and Player Valuations

As part of the initial plots, we have the flows for each of the strategies in the problem. The equations above will adjust for different models in which each player can have any number of strategies. The game “morra” has quite a lot for example.

Table[Plot3D[streamlines[s1,s2][[k]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},PlotRange->Full,ImageSize->Large,PlotLabel->StringJoin["Covariant Flow V: ",strategyString[[k]],"((Subscript[s, 1]),(Subscript[s, 2]))"]],{k,numActiveStrategies}]

As anticipated, the streamlines for the inactive strategies are no longer constants, nor is the streamline for the energy density. We see two-dimensional standing wave structures.

Plot3D[streamlines[s1,s2][[-3]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},ImageSize->Large,PlotLabel->"Energy Density μ((Subscript[s, 1]),(Subscript[s, 2]))",PlotRange->Full]
Plot3D[streamlines[s1,s2][[-2]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])((Subscript[s, 1]),(Subscript[s, 2]))",PlotRange->Full]
Plot3D[streamlines[s1,s2][[-1]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 2])((Subscript[s, 1]),(Subscript[s, 2]))",PlotRange->Full]
Plot3D[streamlines[s1,s2][[-2]]+streamlines[s1,s2][[-1]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])((Subscript[s, 1]),(Subscript[s, 2]))+(Subscript[e, 2])((Subscript[s, 1]),(Subscript[s, 2]))",PlotRange->Full]
Plot3D[streamlines[s1,s2][[-2]]-streamlines[s1,s2][[-1]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])((Subscript[s, 1]),(Subscript[s, 2]))-(Subscript[e, 2])((Subscript[s, 1]),(Subscript[s, 2]))",PlotRange->Full]

Next Iteration and the Covariant Energy as a Test of Convergence

These results are provisional since the PDE is non-linear. To get the next iteration we set the following substitution rule and increment the iteration number by one:

sub=Rule@@@Transpose[{Join[ownershipVariables,preferenceList,{μ[s1,s2]},inactiveOwners],streamlines[s1,s2]}][[All,;;2]]
numberIterations=numberIterations+1

Based on the results so far, we would like to compare the current value with the older value and see if the difference is starting to converge.

ψ1[s1_,s2_]=((-((uI1[s1,s2])^2)-owner[s1,s2].one.owner[s1,s2])/.sub)Exp[-2Φ[preferenceList/.sub]];
ψ2[s1_,s2_]=((-((uI2[s1,s2])^2)-owner[s1,s2].two.owner[s1,s2])/.sub)Exp[-2Φ[preferenceList/.sub]];
v0New[s1_,s2_]=(√(1-ψ1[s1,s2]-ψ2[s1,s2]))(Exp[ν[preferenceList/.sub]]);
deltav0[s1_,s2_]=v0New[s1,s2]-v0Old[s1,s2];

Plot3D[v0Old[s1,s2],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},PlotRange->Full,ImageSize->Large,PlotLabel->"Covariant Flow (Subscript[V, ο])((Subscript[s, 1]),(Subscript[s, 2]))"]
StringJoin[ToString[numberIterations], " iterations"]
Plot3D[deltav0[s1,s2],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},PlotRange->Full,ImageSize->Large,PlotLabel->"(Subscript[δV, ο])((Subscript[s, 1]),(Subscript[s, 2]))"]

If we see this difference go to zero, we have some hope that we have convergence and a consistent overall solution. At this point, if you want to do an iteration, run this section again.Our experience so far with the “known” source choice with no independent variables in the source, is that you may have to run the iterations a number of times. 

We are on iteration number

numberIterations

It is instructive to follow the iterations to get a sense of the structures that emerge. For example, sharp behaviors may indicate that a wave structure with small oscillations is missing and about to be generated. This is especially possible if the boundary condition sets up standing waves that meet inside the region, reflect off of each other and create a sharp boundary. there should be reflections off this interface. Looking at the different iterations is similar to watching a picture in solution slowly develop.

It is also useful to watch the iterations since you can almost anticipate when the solution errors have started to dominate and the solution will start to diverge.

Behavior Along the Streamline s=(1/2)((Subscript[s, 1])+(Subscript[s, 2]))

We make the connection with the harmonic analysis of the previous section by looking at the streamlines along the diagonal (Subscript[s, 1])=(Subscript[s, 2])=s. We start with the energy density and inactive flows. Note that when the inactive flows go to zero, we expect that the energy density is also constant. So to first approximation we might expect that these terms will be constants, which was assumed in the harmonic analysis case.

Plot[streamlines[s,s][[-3]],{s,0,size},AxesLabel->{s,μ},ImageSize->Large,PlotLabel->"μ(s,s)",PlotRange->Full,PlotStyle->Purple]
Plot[streamlines[s,s][[-2]],{s,0,size},AxesLabel->{s,"(Subscript[e, 1])"},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])(s,s)",PlotRange->Full,PlotStyle->Blue]
Plot[streamlines[s,s][[-1]],{s,0,size},AxesLabel->{s,"(Subscript[e, 2])"},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[e, 2])(s,s)",PlotRange->Full,PlotStyle-> Red]

Plot[{streamlines[s,s][[-2]]-streamlines[0,0][[-2]],streamlines[s,s][[-1]]-streamlines[0,0][[-1]]},{s,0,size},AxesLabel->{s,"(Subscript[δe, k])"},ImageSize->Large,PlotLabel->"Covariant Inactive Flow (Subscript[δe, k])(s,s)",PlotRange->Full,PlotStyle-> {Blue,Red},PlotRange->Full,PlotLegends->{"Player 1","Player 2"}]
ParametricPlot[{streamlines[s,s][[-2]],streamlines[s,s][[-1]]},{s,0,size},AxesLabel->{"(Subscript[e, 1])","(Subscript[e, 2])"},PlotStyle->{Black,Thick},PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])(s,s) vs (Subscript[e, 2])(s,s)",ImageSize->Large]

By contrast, the active flows should be harmonic, albeit the harmonics appear to be shifted between players.

Plot[Evaluate[Join[Table[streamlines[s,s][[k]],{k,numActiveStrategies}],{v0[s,s]}]],{s,0,size},PlotLabel->"CoVariant Flow (Subscript[V, a])(s,s)",AxesLabel->{s,"(Subscript[V, a])"},ImageSize->Large,PlotRange->Full,PlotStyle->allStrategyColors,PlotLegends->allStrategyString]

spaceMetric=Table[metric[[i]][[j]],{i,numActiveStrategies},{j,numActiveStrategies}]
contraStreamlines[s1_,s2_]=spaceMetric.Take[streamlines[s1,s2],numActiveStrategies]

Plot[Evaluate[Join[Table[contraStreamlines[s,s][[k]],{k,numActiveStrategies}],{v0[s,s]}]],{s,0,size},PlotLabel->"ContraVariant Flow (V^a)(s,s)",AxesLabel->{s,"(V^a)"},ImageSize->Large,PlotRange->Full,PlotStyle->allStrategyColors,PlotLegends->allStrategyString]

For comparison, the solution from single ownership model:

Plot[Evaluate[Join[Table[-flowSolution[s][[k]]+contraStreamlines[s,s][[k]],{k,numActiveStrategies}],{-flowSolution[s][[-1]]+v0[s,s]}]],{s,0,sMax},PlotLabel->"ContraVariant Flow comparison ((δV)^a)(s,s)",AxesLabel->Automatic,PlotLegends->Join[strategyString2,strategyString1,{"time"}],PlotStyle->Join[colors2,colors1,{Purple}],ImageSize->Large,PlotRange->All]

Now back to two owners:

Manipulate[Plot[Evaluate[Join[Table[contraStreamlines[s+shift,s-shift][[k]],{k,numActiveStrategies}],{v0[s+shift,s-shift]}]],{s,Min[Abs[shift],size-Abs[shift]],Max[Abs[shift],size-Abs[shift]]},PlotLabel->"ContraVariant Flow (V^a)(s+shift,s-shift)",AxesLabel->{"s","(V^a)"},PlotRange->Full,PlotStyle->allStrategyColors,PlotLegends->allStrategyString],{{shift,0},-(size/2),(size/2)},SaveDefinitions->True]

Manipulate[TableForm[{{Plot[{streamlines[s+shift,s-shift][[-2]]-streamlines[0,0][[-2]],streamlines[s,s+shift][[-1]]-streamlines[0,0][[-1]]},{s,Min[Abs[shift],size-Abs[shift]],Max[Abs[shift],size-Abs[shift]]},AxesLabel->{"s","(Subscript[δe, k])"},PlotLabel->"Covariant Inactive Flow (Subscript[δe, k])(s+shift,s-shift)",PlotRange->Full,PlotStyle-> {Blue,Red},PlotRange->Full,PlotLegends->{"Player 1","Player 2"}],Plot[streamlines[s+shift,s-shift][[-3]],{s,Min[Abs[shift],size-Abs[shift]],Max[Abs[shift],size-Abs[shift]]},AxesLabel->{"s",μ},PlotLabel->"μ(s+shift,s-shift)",PlotRange->Full,PlotStyle->Purple,PlotRange->Full]},{
Plot[streamlines[s+shift,s-shift][[-2]],{s,Min[Abs[shift],size-Abs[shift]],Max[Abs[shift],size-Abs[shift]]},AxesLabel->{"s","(Subscript[e, 1])"},PlotLabel->"Covariant Inactive Flow (Subscript[e, 1])(s+shift,s-shift)",PlotRange->Full,PlotStyle->Blue,PlotRange->Full],
Plot[streamlines[s+shift,s-shift][[-1]],{s,Min[Abs[shift],size-Abs[shift]],Max[Abs[shift],size-Abs[shift]]},AxesLabel->{"s","(Subscript[e, 2])"},PlotLabel->"Covariant Inactive Flow (Subscript[e, 2])(s+shift,s-shift)",PlotRange->Full,PlotStyle-> Red,PlotRange->Full]}}],{{shift,0},-(size/2),(size/2)},SaveDefinitions->True]

Preferences and Strategic Efforts

streamlines[size,size][[2numActiveStrategies+1]]
streamlines[size,0][[2numActiveStrategies+1]]
xSolution[size][[numActiveStrategies+1]]
size flowSolution[0][[numActiveStrategies+1]]

Plot[Evaluate[Table[streamlines[s,s][[k]],{k,numActiveStrategies+1,2numActiveStrategies+1}]],{s,0,size},PlotLabel->"Preference (y^a)(s,s)",AxesLabel->{s,"(y^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]
LogPlot[Evaluate[Table[Exp[streamlines[s,s][[k]]],{k,numActiveStrategies+1,2numActiveStrategies}]],{s,0,size},PlotLabel->"Strategic Effort (E^((y^a) (s, s)))",AxesLabel->{s,"(y^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]

Plot[Evaluate[Table[-xSolution[s][[k-numActiveStrategies]]+streamlines[s,s][[k]],{k,numActiveStrategies+1,2numActiveStrategies+1}]],{s,0,sMax},PlotLabel->"ContraVariant Flow comparison ((δy)^a)(s,s)",AxesLabel->Automatic,PlotLegends->allStrategyString,PlotStyle->allStrategyColors,ImageSize->Large,PlotRange->Full]

Table[Plot3D[preferenceNumerical[s1_,s2_]=streamlines[s1,s2][[k]],{s1,s2}∈Ω1,AxesLabel->{(Subscript[s, 1]),(Subscript[s, 2])},PlotRange->Full,ImageSize->Large,PlotLabel->StringJoin["preference (y^a): ",allStrategyString[[k-numActiveStrategies]],"((Subscript[s, 1]),(Subscript[s, 2]))"]],{k,numActiveStrategies+1,2numActiveStrategies+1}]

Forced Harmonic Disruptive Influences

Introduction to harmonic disruptions

We deal here with the possibility that there may be initial impulses on the system; we want to study the behavior of the system due to these impulses. We follow the technique from other areas, including electrical engineering, to look at the system in terms of its phasor components (Fourier components if you like): forced harmonics. A study of these forced harmonics over all frequencies is equivalent to a study of what happens to a system under impulses. What will be of interest is whether there are resonant structures, so that a hit to the system will get amplified. These are observed.

This section is not complete because we don’t yet have a way to Solve the linear Laplace equation in more than 3 dimensions. We present some illustrative examples however in 2 and 3 dimensions to get a feel for the issues. In the previous sections we dealt with two ways to generate harmonic behaviors: payoff contributions and seasonal payoff contributions. One additional type of stationary flow results from perturbing the system with an impulse at some initial point of time. These flows in general will be in a frame in which the metric is time dependent. In general, we show that an initial impulse or perturbation can be expressed as a linear combination over (an infinite) set of harmonic contributions. We study this full set of harmonic perturbations on the system.

The way we approach harmonic disruptive influences is to start with our initial active metric (g^(ab)) and inactive metric (Subscript[γ, jk]) structures in the holonomic stationary frame, in which these metric components are independent of time. Let the coordinates for this frame be (y^a) . 

Realistic behaviors will have metric components that are time dependent. Assume that we have such a metric and choose the covariant gauge with coordinates (x^a). It follows that these coordinates satisfy the wave equation. We write this wave equation in the holonomic stationary frame

					({{(x^(Overscript[υ, ¯]))=(Subscript[(x^(Overscript[υ, ¯])), ω])(e^(iωt))}, {(ω^2)(Subscript[(x^(Overscript[υ, ¯])), ω])−(Subscript[(Overscript[g, ¯]), οο])((Overscript[g, ¯])^(υυ'))(∂_(υ'))(∂_υ)(Subscript[(x^(Overscript[υ, ¯])), ω])=0}, {(δ^(υυ'))(∂_(υ'))(∂_υ)(Subscript[(x^(Overscript[υ, ¯])), ω])+(ω^2)(e^(−2(ν−Φ)))(Subscript[(x^(Overscript[υ, ¯])), ω])=0}, {(δ^(υυ'))(∂_(υ'))(∂_υ)(Subscript[(x^(Overscript[υ, ¯])), ω])+(ω^2)(e^(−2((n−1)/(n−2))ν))(Subscript[(x^(Overscript[υ, ¯])), ω])=0}})

Note that we can make the frequency for each coordinate different, corresponding to different impulses along each strategic preference direction. These are Laplace equations that we can solve numerically. The challenge may be that Mathematica may not have a general purpose solver using Dirichlet conditions for dimensions higher than three. 

There are some initializations that will make our life somewhat easier. To deal with the “current models” below, these initializations have been recorded in an earlier section under initializations “Useful Methods” at the start of the notebook on student projects.

Basic Laplace Equation in 2 and 3 dimensions

We start by restricting attention to 2 dimensions, and redefine some of the constants

ySize=11
yScale=1
νScale=1+1/3
ωP=2
scaling=(1/10)
frameModel[(((Subscript[y, 1])-(Subscript[y, 2]))/yScale)]
ContourPlot[Exp[-νScale frameModel[(((Subscript[y, 1])-(Subscript[y, 2]))/yScale)]],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},PlotLegends->Automatic]
Plot3D[Exp[-νScale frameModel[(((Subscript[y, 1])-(Subscript[y, 2]))/yScale)]],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},PlotRange->Full,AxesLabel->Automatic]

We only have two strategy directions and we can use Mathematica to solve the Laplace equation

prefyS={(Subscript[y, 1]),(Subscript[y, 2])}
Ω2=Cuboid[{0,0},{ySize,ySize}]
xNu=NDSolveValue[{harmonic[prefyS,ωP Exp[-νScale frameModel[(((Subscript[y, 1])-(Subscript[y, 2]))/yScale)]],x],DirichletCondition[x@@prefyS==scaling (Subscript[y, 1]),True]},x,prefyS∈Ω2,MaxSteps->∞]
Plot3D[xNu[(Subscript[y, 1]),(Subscript[y, 2])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize}]

Note that the solution is well behaved.

3 Dimension example

We can also solve the problem in 3 dimensions. We imagine there three strategies, the first two owned by player 1 and the third by player 2:

prefyS={(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])}
Ω2=Cuboid[{0,0,0},{ySize,ySize,ySize}]
xNu=NDSolveValue[{harmonic[prefyS,ωP Exp[-νScale frameModel[(((((Subscript[y, 1])+(Subscript[y, 2]))/2)-(Subscript[y, 3]))/yScale)]],x],DirichletCondition[x@@prefyS==scaling (Subscript[y, 2]),True]},x,prefyS∈Ω2,MaxSteps->∞]
Plot3D[xNu[(Subscript[y, 1]),(ySize/2),(Subscript[y, 3])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 3]),0,ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-3"]
Plot3D[xNu[(Subscript[y, 1]),(Subscript[y, 2]),(ySize/2)],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-2"]

Again the solution is well behaved over the full range. We do see the effect of the drift potential in distorting the harmonics. We have used the Dirichlet condition where we imagine the solution x is for the (Subscript[y, 2]). The value of x on the boundary is (Subscript[y, 2]) , which we scaled for more generality. 

We illustrate the effects in 3 dimensional density and slice density plots

DensityPlot3D[xNu[(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},{(Subscript[y, 3]),0,ySize},ColorFunction->Hue,OpacityFunction->Automatic,PlotLegends->Automatic,AxesLabel->Automatic]

SliceDensityPlot3D[xNu[x,y,z],{x-y==0 ,x+y==ySize ,z==(ySize/2)},{x,0,ySize},{y,0,ySize},{z,0,ySize},ColorFunction->Hue,PlotLegends->Automatic,AxesLabel->Automatic]

We have applied the Dirichlet Boundary condition everywhere, which enforces the stable behavior.

3 Dimensional Example with Dirichlet only on the bottom and sides for player 2

We take a first step in moving towards the method of lines solution by not imposing the Dirichlet condition when (Subscript[y, 2])==0.

xNuMOL=NDSolveValue[{harmonic[prefyS,ωP Exp[-νScale frameModel[(((((Subscript[y, 1])+(Subscript[y, 2]))/2)-(Subscript[y, 3]))/yScale)]],x],DirichletCondition[x@@prefyS==scaling (Subscript[y, 2]),(Subscript[y, 2])<ySize]},x,prefyS∈Ω2,MaxSteps->∞]
Plot3D[xNuMOL[(Subscript[y, 1]),(ySize/2),(Subscript[y, 3])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 3]),0,ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-3"]
Plot3D[xNuMOL[(Subscript[y, 1]),(Subscript[y, 2]),(ySize/2)],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-2"]

Note that the solution along (Subscript[y, 2]) is different from before. It is not constrained by the top boundary. This changes the density plots as well.

DensityPlot3D[xNuMOL[(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),0,ySize},{(Subscript[y, 3]),0,ySize},ColorFunction->Hue,OpacityFunction->Automatic,PlotLegends->Automatic,AxesLabel->Automatic]

SliceDensityPlot3D[xNuMOL[x,y,z],{x-y==0 ,x+y==ySize ,z==(ySize/2)},{x,0,ySize},{y,0,ySize},{z,0,ySize},ColorFunction->Hue,PlotLegends->Automatic,AxesLabel->Automatic]

Method of Lines for the 3 dimensional example

We revisit the above examples using the Method of Lines midway up the value for the coordinate (Subscript[y, 2]).  We allow the possibility to change the parameters.

νScale =1
yScale=1
ySize=11
scaling=(1/10)

In using the Method of Lines, we ran into a problem with Mathematica; we were not successful using the indexed variables (Subscript[y, k]) in NDSolve. We do a change of variables nDVariables, that works and create the region.

nDVariables
region

We write out the equations as we did in Vol 2, where we assume with the Method of Lines we will have a “time” direction (Subscript[y, 2]) with initial value conditions abnd periodic boundary conditions on the other two directions:

equMOL={harmonic[prefyS,ωP Exp[-νScale frameModel[(((((Subscript[y, 1])+(Subscript[y, 2]))/2)-(Subscript[y, 3]))/yScale)]],w],(w@@prefyS/.(Subscript[y, 2])->(ySize/2))==(ySize/2) scaling Cos[((2Pi (Subscript[y, 1]))/ySize)]Cos[((2Pi (Subscript[y, 3]))/ySize)],(D[w@@prefyS,(Subscript[y, 2])]/.(Subscript[y, 2])->(ySize/2))==scaling ,w[0,(Subscript[y, 2]),(Subscript[y, 3])]==w[ySize,(Subscript[y, 2]),(Subscript[y, 3])],w[(Subscript[y, 1]),(Subscript[y, 2]),0]==w[(Subscript[y, 1]),(Subscript[y, 2]),ySize]}/.nDVariables

We use the same parameters for NDSolve that we used in Vol. 2

xNuMOL2=NDSolveValue[equMOL,w,Sequence@@region,MaxSteps->∞,Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MinPoints"->21,"DifferenceOrder"->"Pseudospectral"}}]
Plot3D[xNuMOL2[(Subscript[y, 1]),(ySize/2),(Subscript[y, 3])],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 3]),0,ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-3"]
Plot3D[xNuMOL2[(Subscript[y, 1]),(Subscript[y, 2]),(ySize/2)],{(Subscript[y, 1]),0,ySize},{(Subscript[y, 2]),((1/2)-.1)ySize,((1/2)+.1)ySize},AxesLabel->Automatic,PlotRange->Full,PlotLabel->"1-2"]

We see the downside of the Method of Lines: The solution has a small region in which it gives plausible results since outside this region the solutions become arbitrarily large. We also took the starting point in the middle of the region and have to assume something about the periodic nature in that region. For these reasons we look for some better method. The upside of the method is that it is general purpose and gives us some idea of what is going on.

Special Case Spectral Decomposition: frame effects depend only on player effort (Subscript[r, k])

The wave equation can be decomposed in a spectral decomposition, which helps in special cases. We consider one such case in which the frame effects depend only on the player efforts (Subscript[r, J]). We make a transformation to a coordinate frame to reflect this:

				({{(Subscript[r, J])=(1/(Subscript[n, J]))(∑_(j[J]))(Subscript[y, j[J]])}, {(Subscript[q, j[J]])=(Subscript[y, j[J]])−(Subscript[r, J])}, {(∑_(j[J]))(Subscript[q, j[J]])=0⇒(Subscript[q, (Subscript[n, J])])≡−(∑_(j[J]≠(Subscript[n, J])))(Subscript[q, j[J]])}})

We could write an orthogonal transformation, but for our purposes, this won’t change the final result. We assume that the frame effects depend only on the efforts, so that the spectral decomposition in the relative variables φ=(Subscript[(Overscript[φ, ˜]), (Subscript[j, 0])])(e^(λ(Subscript[q, (Subscript[j, 0])]))) generates  only the solutions with λ=0. We have the following wave equations:

				({{(∑_K)(1/(Subscript[n, K]))(((∂^2)(Subscript[(Overscript[φ, ˜]), j[J]]))/(∂((Subscript[r, K])^2)))+ω((((Subscript[r, 1]),...,(Subscript[r, N])))^2)(Subscript[(Overscript[φ, ˜]), j[J]])=0}, {(∑_K)(1/(Subscript[n, K]))(((∂^2)(Subscript[(Overscript[φ, ˜]), 0]))/(∂((Subscript[r, K])^2)))+ω((((Subscript[r, 1]),...,(Subscript[r, N])))^2)(Subscript[(Overscript[φ, ˜]), 0])=0}})

These are the same equation and differ only in the boundary conditions:

				({{(Subscript[(Overscript[φ, ˜]), j[J]])|=1⇒(Subscript[(Overscript[φ, ˜]), j[J]])≡(Subscript[(Overscript[φ, ˜]), 0])}, {(Subscript[φ, j[J]])|=(Subscript[q, j[J]])}, {φ|=(Subscript[(Overscript[φ, ˜]), (Subscript[r, J])])|=(Subscript[r, J])}})

Based on these equations, we recover the required coordinate scalars:

				(x^(a[J]))=(Subscript[(Overscript[φ, ˜]), (Subscript[r, J])])+(Subscript[(Overscript[φ, ˜]), 0])(Subscript[q, a[J]])

These equations will work for all of our models, require only that the frame effects depend on the efforts. For two players we have three wave equations to solve based on their different wave equations.

Finally, note we can temporarily scale our efforts so that we get a true Laplacian

				(Subscript[r, J])≡((Subscript[(Overscript[r, ¯]), J])/(Sqrt[(Subscript[n, J])]))
So we start by defining the variables and ranges

prefyS=prefy
r1=((Tr[prefy.one])/(Tr[one]))
r2=((Tr[prefy.two])/(Tr[two]))
qOne=Simplify[one.prefy-r1 Table[1,{k,numActiveStrategies+1}].one]
qTwo=Simplify[two.prefy-r2 Table[1,{k,numActiveStrategies+1}].two]

rbar1=((Tr[prefy.one])/(Sqrt[Tr[one]]))
rbar2=((Tr[prefy.two])/(Sqrt[Tr[two]]))

prefyR=(qOne+((Subscript[r, 1])/(Sqrt[Tr[one]]))).one+(qTwo+((Subscript[r, 2])/(Sqrt[Tr[two]]))).two

Simplify[ν[prefyR]]
Simplify[Φ[prefyR]]
Simplify[Φ[prefyR]-ν[prefyR]]
prefySR={(Subscript[r, 1]),(Subscript[r, 2])}
ωP0=setForcedHarmonicFrequency
xScale=1
ySize=setMaxPreferenceForHarmonicCalculation
ωP=ωP0  Exp[Simplify[Φ[prefyR]-ν[prefyR]]]
harmonic[prefySR,ωP,w]

So now we compute the relative strategy solutions

ySize1=(Sqrt[Tr[one]])ySize
ySizeLow1=- (Sqrt[Tr[one]])ySize
ySize2=(Sqrt[Tr[two]])ySize
ySizeLow2=-(Sqrt[Tr[two]])ySize
Ω2=Cuboid[{ySizeLow1,ySizeLow2},{ySize1,ySize2}]//N
rNu0=NDSolveValue[{harmonic[prefySR,ωP,w],DirichletCondition[w@@prefySR==xScale,True]},w,prefySR∈Ω2,MaxSteps->∞]
Plot3D[rNu0[(Subscript[r, 1]),(Subscript[r, 2])],{(Subscript[r, 1]),ySizeLow1,ySize1},{(Subscript[r, 2]),ySizeLow2,ySize2},AxesLabel->Automatic,PlotLabel->"(Subscript[ϕ, 0])"]

We change the boundary condition for player 1 solutions

rNu1=NDSolveValue[{harmonic[prefySR,ωP,w],DirichletCondition[w@@prefySR==xScale (Subscript[r, 1]),True]},w,prefySR∈Ω2,MaxSteps->∞]
Plot3D[rNu1[(Subscript[r, 1]),(Subscript[r, 2])],{(Subscript[r, 1]),ySizeLow1,ySize1},{(Subscript[r, 2]),ySizeLow2,ySize2},AxesLabel->Automatic,PlotLabel->"(Subscript[ϕ, 1])"]

For Player 2:

rNu2=NDSolveValue[{harmonic[prefySR,ωP,w],DirichletCondition[w@@prefySR==xScale (Subscript[r, 2]),True]},w,prefySR∈Ω2,MaxSteps->∞]
Plot3D[rNu2[(Subscript[r, 1]),(Subscript[r, 2])],{(Subscript[r, 1]),ySizeLow1,ySize1},{(Subscript[r, 2]),ySizeLow2,ySize2},AxesLabel->Automatic,PlotLabel->"(Subscript[ϕ, 2])"]

These are combined to form the new coordinates, defining first the numerical preferences from the previous section

preferenceNumerical[s1_,s2_]=Table[streamlines[s1,s2][[k+numActiveStrategies]],{k,numActiveStrategies+1}];

The results are:

xNum[pref_]:=((qOne rNu0[rbar1,rbar2]+((rNu1[rbar1,rbar2])/(Sqrt[Tr[one]]))).one+(qTwo rNu0[rbar1,rbar2]+((rNu2[rbar1,rbar2])/(Sqrt[Tr[two]]))).two+(metric+one+two).Table[1,{k,numActiveStrategies+1}]prefy[[-1]])/.Rule@@@Transpose[{prefy,Table[pref[[k]],{k,numActiveStrategies+1}]}]

xNumPref[s1_,s2_]=xNum[preferenceNumerical[s1,s2][[1]],preferenceNumerical[s1,s2][[2]],preferenceNumerical[s1,s2][[3]],preferenceNumerical[s1,s2][[4]],preferenceNumerical[s1,s2][[5]]];
xNumPref[s1_,s2_]=xNum[preferenceNumerical[s1,s2]];

Table[Plot3D[preferenceNumerical[s1,s2][[k]],{s1,0,sMax},{s2,0,sMax},AxesLabel->Automatic,PlotRange->Full,PlotLabel->StringJoin["(y^a) a=",allStrategyString[[k]]]],{k,numActiveStrategies+1}]
Table[Plot3D[xNumPref[s1,s2][[k]],{s1,0,sMax},{s2,0,sMax},AxesLabel->Automatic,PlotRange->Full,PlotLabel->StringJoin["(x^a) a=",allStrategyString[[k]]]],{k,numActiveStrategies+1}]

Plot[Evaluate[Table[streamlines[s,s][[k]],{k,numActiveStrategies+1,2numActiveStrategies+1}]],{s,0,size},PlotLabel->"Preference (y^a)(s,s)",AxesLabel->{s,"(y^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]
Plot[Evaluate[Table[xNumPref[s,s][[k]],{k,numActiveStrategies+1}]],{s,0,size},PlotLabel->"Preference (x^a)(s,s)",AxesLabel->{s,"(x^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]
LogPlot[Evaluate[Table[Exp[xNumPref[s,s][[k]]],{k,numActiveStrategies}]],{s,0,size},PlotLabel->"Strategic Effort (E^((x^a) (s, s)))",AxesLabel->{s,"(x^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]

With Time Dependence Included

2.Pi/ωP0
compNumPrefSin[s_]=Table[Sin[ωP0 streamlines[s,s][[1+2numActiveStrategies]]]xNumPref[s,s][[k]],{k,numActiveStrategies}];
compNumPrefCos[s_]=Table[Cos[ωP0 streamlines[s,s][[1+2numActiveStrategies]]]xNumPref[s,s][[k]],{k,numActiveStrategies}];
Plot[Evaluate[compNumPrefSin[s]],{s,0,size},PlotLabel->"Preference sin(ωt)(x^a)(s,s)",AxesLabel->{s,"sin(ωt)(x^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]
Plot[Evaluate[compNumPrefCos[s]],{s,0,size},PlotLabel->"Preference cos(ωt)(x^a)(s,s)",AxesLabel->{s,"cos(ωt)(x^a)"},ImageSize->Large,PlotStyle->allStrategyColors,PlotLegends->allStrategyString,PlotRange->Full]

Eigensystem

The system depends on the eigenvalues and eigenfrequencies of the Laplacian, which are given below. Just as in physical systems, when looking at the frequency decomposition, the new preferences as a function of frequency may have large amplitudes near resonances of the system. The eigensystems below are just such resonances.

((ωP0 )^2)//N
Plot3D[-rNu0[(Subscript[r, 1]),(Subscript[r, 2])],{(Subscript[r, 1]),ySizeLow1,ySize1},{(Subscript[r, 2]),ySizeLow2,ySize2},AxesLabel->Automatic,PlotLabel->"-(Subscript[ϕ, 0])"]

{prefVals,prefFuns}=NDEigensystem[{-(( Exp[Simplify[-Φ[prefyR]+ν[prefyR]]])/.Rule@@@Transpose[{{(Subscript[r, 1]),(Subscript[r, 2])},{x,y}}])Laplacian[w[x,y],{x,y}],DirichletCondition[w[x,y]==0,True]},w[x,y],{x,y}∈Ω2,20];

prefVals
Table[Plot3D[prefFuns[[i]],{x,y}∈Ω2,PlotRange->All,PlotLabel->prefVals[[i]],PlotTheme->"Minimal",AxesLabel->{"(Subscript[r, 1])","(Subscript[r, 2])"}],{i,Length[prefVals]}]

Conclusions

Lessons Learned

I should not be shy but call the course “A computational approach to decision process theory”. This is where I got to at the end. This is a new approach to a  theory of decision processes.

The focus on games and their harmonic behaviors was the right thing to do and moreover, it was crucial to provide the "initializations" that everyone could use in their own notebooks.

If I taught this course again,

I would not introduce the symmetric form of game theory but go directly to the decision process theory form: having two alternative forms was an unnecessary confusion.

I would remove from the list of initializations functions those I don't intend the students to use such as things referring to conjugate, or dropping the time part of the payoff metric and calling the result "the payoff matrix". Again this caused unnecessary confusion.

I would drill studentd more on what I currently call the player valuations, and not use the concept hedge strategy (again, since it is confusing). The player valuation, along with the payoff fields, form two of the key dynamic mechanisms (yes there are others in the theory as well)

I would do a better job introducing the concept of The Nash equilibrium, especially with regard to games such as the prisoner's dilemma and the various games of chicken. What the theory does that is new is, for games that have two or more Nash equilibria, provides a systematic way of combining them

I was caught off guard on other possible "harmonic" behaviors for a general antisymmetric matrix and their relationship to Nash equilibria.

There are harmonics, which correspond to imaginary eigenvalues of the mixed tensor and oscillations around a stable equilibrium point.

A special case of harmonics is when the null space consists of a single timelike vector with zero space components: this is the co-moving equilibrium point. The preferences don't grow in this case but stay bounded.

However, there was a case that escaped my attention, corresponding to an unstable equilibrium point. In this case there may not be a null space, or only a null space with spacelike vectors. There will be real eigenvalues so the behavior might be like the hyperbolic behaviors in classical physics. There seems to be nothing wrong with such behaviors. I just note that starting with a timelike flow vector, the length of the flow will be maintained by the equations.

These three behaviors should be introduced early in the course and then supported with homework exercises.
Changes focus on null space of the payoff and the disruptive startup of the flow; the null space and eigenvalues give clues for the type of fixed points we can have.

I believe the exercises helped the students grasp the concepts

Each student presented their final project to the whole class. This was particularly instructive as it gave students a view of how others had learned and so everyone go pulled up by their collective bootstraps.

The game of chicken was chosen by several students and was particularly interesting and relevant. It showed that the null space need not be a Nash equilibrium. This pulled together many of the components of game theory and harmonics.

The harmonics provided an essential insight into decision process theory and show the deviations from Bayesian Game Theory.

For students whose only exposure to computing was Java and similar languages, they found that Mathematica was somewhat brittle and unhelpful when programming mistakes were made, so this caused some grousing. However, most of the issues were quite trivial and easy to fix once identified.

Because of the length of the quarter system, I didn’t have time to do more than just sketch the possibilities of player ownership and frame effects. These are two additional dynamic mechanims that I think might have been of interest if there had  been more time.