# My model has not finished evaluating in more than a day so can't test if it works, what is wrong with it?

Hello I have this modified monocentric city model that I have been working on for a while.

Sorry for the long post, I wanted to add the whole code, because it's all connected.

My suspicion is that it has something to do with the p (price) and q (demand) functions [economics]. My previous version used a simpler price and demand functions, but since I changed them, I have not been able to get it to work. I am no mathematical genius, but I have tried so many things that I'm about to give up hope on this model. I hope one of you can help me figure out the problem.

(* Cost share of housing capital in housing production function *)
\[Beta] = 0.6;
(* Scaling on housing production function *)
g = 0.0005;
(* Radians available for construction *)
(* benchmark \[Theta] is 3 *)
\[Theta] = 2*Pi;
(* elasticity of income *)
\[Gamma] = 0.2;
(* price elasticity of housing *)
\[Epsilon] = -0.6;
(* constant in demand function *)
\[Omega] = 1;

p[x_, y_, t_, ta_, f_, fa_,
u_] := (E^\[Omega]*(\[Epsilon] +
1)*((y - (f + fa) - (t + ta)*x)^(1 - \[Gamma])/(1 - \[Gamma]) -
u))^(1/(\[Epsilon] + 1));
q[x_, y_, t_, ta_, f_, fa_,
u_] := (1 - \[Gamma])^(\[Gamma]/(1 - \[Gamma]))*E^\[Omega]*
p[x, y, t, ta, f, fa,
u]^(\[Epsilon] +
1)*((E^\[Omega]*
p[x, y, t, ta, f, fa, u]^(\[Epsilon] + 1))/(\[Epsilon] + 1) +
u)^(\[Gamma]/(1 - \[Gamma]));
S[x_, y_, t_, ta_, f_, fa_,
u_] := (1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] -
1));
r[x_, y_, t_, ta_, f_, fa_,
u_] := (p[x, y, t, ta, f, fa,
u] (g)) ((1/(p[x, y, t, ta, f, fa,
u] (\[Beta]) (g)))^(\[Beta]/(\[Beta] - 1))) -
1 ((1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] - 1)));
h[x_, y_, t_, ta_, f_, fa_,
u_] := (g) S[x, y, t, ta, f, fa, u]^(\[Beta]);
Density[x_, y_, t_, ta_, f_, fa_, u_] :=
h[x, y, t, ta, f, fa, u]/q[x, y, t, ta, f, fa, u];
L[xbar_, y_, t_, ta_, f_, fa_, u_] := \[Theta]*\!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$xbar$$]$$x*\ Density[x, y, t, ta, f, fa, u] \[DifferentialD]x$$\);
pbar[xbar_, y_, t_, ta_, f_, fa_, u_,
pop_] := (\[Theta]/pop)*(1/30000)*\!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$xbar$$]$$x* Density[x, y, t, ta, f, fa, u]* p[x, y, t, ta, f, fa, u] \[DifferentialD]x$$\);
qbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$xbar$$]$$x* Density[x, y, t, ta, f, fa, u]* q[x, y, t, ta, f, fa, u] \[DifferentialD]x$$\);
hbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$xbar$$]$$x* Density[x, y, t, ta, f, fa, u]* h[x, y, t, ta, f, fa, u] \[DifferentialD]x$$\);
xhat[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$xbar$$]$$x^2* Density[x, y, t, ta, f, fa, u] \[DifferentialD]x$$\);
rev[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] :=
pop*(fa + ta*xhat[xbar, y, t, ta, f, fa, u, pop]);
(*fa[xbar_,y_,t_,ta_,f_,fa_,u_,pop_] := \
((rev[xbar,y,t,ta,f,fa,u,pop]/pop)/(ta*xhat[xbar,y,t,ta,f,fa,u,pop]));
xhatnew[xbar_,y_,t_,ta_,f_,fa_,u_,pop_]:= \
((rev[xbar,y,t,ta,f,fa,u,pop]/pop)-fa[xbar,y,t,ta,f,fa,u,pop])/ta;*)

Population = 800000;
y = 70000;
t = 350;
f = 3000;
ta = 250;
fa = 1000;
ra = 45000;

Lu[u_] := (solution =
FindRoot[{r[xbar, y, t, ta, f, fa, u] - ra == 0}, {xbar, 10},
AccuracyGoal -> 1]; xbar /. solution);

fu[u_] := N[L[Evaluate[Lu[u]], y, t, ta, f, fa, u] - Population];

findRoot[fun_, lower_, upper_, eps_] :=
Module[{ul = lower, uh = upper, count = 0},
While[And[Abs[fun[ul]] > eps, count < 10],
If[fun[(uh + ul)/2] > 0, ul = (uh + ul)/2, uh = (uh + ul)/2];
count = count + 1]; Print["Iterations: ", count,
" Population: ", Population, " people",
" utility: ", N[ul],
" City boundary (xbar): ", Lu[N[ul]], " km",
" Average house price (pbar): ",
pbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros",
" Price at 5 km from the centre p(5): ",
p[5, y, t, ta, f, fa, N[ul]]/30000, " euros",
" Price at 10 km from the centre p(10): ",
p[10, y, t, ta, f, fa, N[ul]]/30000, " euros",
" Average density: ", 2*Population/(\[Theta]*Lu[N[ul]]^2),
" people per ?",
" Average floor to area ratio hbar: ",
hbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population],
" Density at 5 km from the centre Density(5): ",
Density[5, y, t, ta, f, fa, N[ul]], " people per ?",
" Density at 20 km from the centre Density(20): ",
Density[20, y, t, ta, f, fa, N[ul]], " people per ?",
" Average distance to the centre xhat: ",
xhat[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " km",
" Tax revenue from transport: ",
rev[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros",
Plot[{Density[x, y, t, ta, f, fa, N[ul]]}, {x, 0, Lu[N[ul]]}]]];
findRoot[fu, 3000, 10000, 1];


It just keeps running. I don't know how long, cause after a day I stopped it.

Extra details:

The p is derived from this indirect utility function (log utility as opposed to my previous Cobb-Douglas utility) through simple equation solving assuming that u is constant:

u = (-e^(\[Omega]*∙))*(p^(\[Epsilon] + 1)/(\[Epsilon] + 1)) +
y^(-\[Gamma] + 1)/(-\[Gamma] + 1)

y is to be replaced by y-(fa+f)-(ta+t)x


The q is derived by taking the u functionand putting y-(fa+f)-(ta+t)x on the left side:

y - (f + fa) - (t + ta)*
x = ((-\[Gamma] + 1)*(u +
e^(\[Omega]*∙)*(p[x, y, t, ta, f, fa,
u]^(\[Epsilon] + 1)/(\[Epsilon] + 1))))^(1/(-\[Gamma] + 1))


and differentiating with respect to p, which should be the equation above for q.

Any help would be greatly appreciated. Also help about the process of troubleshooting this.

Here a couple of changes/suggestions.

Using a decimal dot for all constants forces most algorithms to switch from symbolic to numeric routines with machine floating point arithmetic; that's typically good for performance:

Population = 800000.;
y = 70000.;
t = 350.;
f = 3000.;
ta = 250.;
fa = 1000.;
ra = 45000.;


Localize xbar to bulletproof the function Lu. (This was not an issue here, but xbar appearing in blue rang some alarm bells when I first saw it.)

Lu[u_] := Block[{xbar},
xbar /. FindRoot[{r[xbar, y, t, ta, f, fa, u] == ra}, {xbar, 10}]
];


Use numeric integration instead of symbolic integration:

L[xbar_, y_, t_, ta_, f_, fa_, u_] := \[Theta] NIntegrate[x Density[x, y, t, ta, f, fa, u], {x, 0, xbar}];


Here I removed the Evaluate which got stuck because the symbolic evaluation of L failed. Using the pattern fu[u_?NumericQ] makes sure that nobody accidentally calls fu with symbolic inputs.

ClearAll[fu];
fu[u_?NumericQ] := L[Lu[u], y, t, ta, f, fa, u] - Population;


Now I am able to evaluate the following in reasonable time:

fu[3000.]
fu[10000.]


-793799.

-800000. - 0.934714 I

There will be plenty of errors due to diverging(?) integrals. Also note that 3000. and 10000. were the lower and upper bound for your bisection root finding algorithm. But your algorithm expects that fu[3000.]>0 and that fu[10000.]>0. As you can see, none of that is correct. Find better bounds (if there are any) and start again.

• I can't believe it.. it works, and for me it gives sort of valid results in seconds. I may have to tweak the numbers a bit, but it totally works. I spent MONTHS on this. Thank you so much! May 23 at 15:50
• Yeah, you seemed a bit desperate. I am glad that I could help. God luck with your experiments! May 23 at 15:51