Numerical Solution of System of Difference/Differential Algebraic Equations in Maxima

Question

I begin with a system of difference equations (an oversimplified Solow-Romer economic model):

Where the t subscripts indicate discrete time, such as Y[t=0], Y[t=1], Y[t=2], ...

Specifically:

eq1: Y[t] = A[t]*K[t]^(1/3)*Ly[t]^(2/3);
eq2: K[t+1] - K[t] = s*Y[t] - d*K[t];
eq3: A[t+1] - A[t] = z*A[t]*La[t];
eq4: Ly[t] + La[t] = L;
eq5: La[t] = l*L;

Endogenous variables (unknowns), and their initial conditions:
    Y[t]          Y[0] = 9663.8253
    K[t]          K[0] = 100.0
    A[t]          A[0] = 100.0
    Ly[t]        Ly[0] = 95.0
    La[t]        La[0] = 5.0

Exogenous variables (givens):
    s: 0.15;
    d: 0.07;
    z: 0.02;
    l: 0.05;
    L: 100.0;

This is 5 equations in 5 unknowns. "Solving" the system numerically is trivial in practice: you just start at t=0 with the initial conditions, calculate K[1] and A[1] from the difference equations, and then calculate Y[1] from that.

Despite its trivial nature, I haven't been able to determine how to actually do so and plot the resulting curves in Maxima.

I'm perfectly comfortable with a differential equation approach (really differential-algebraic) if that is more conducive to Maxima's capabilities. That should be equivalent in a numerical solution anyway:

That is:

eq1: Y(t)=A(t)*K(t)^(1/3)*Ly(t)^(2/3);
eq2: diff(K(t),t) = s*Y(t)-d*K(t);
eq3: diff(A(t),t) = z*A(t)*La(t);
eq4: Ly(t)+La(t) = L;
eq5: La(t) = l*L;

But, again, I don't see a way to numerically solve and plot this system with Runge-Kutta or the other built-in solvers (and this is true even though the algebraic equations above can be easily rewritten in the form 0=f(Y,A,K,Ly,La)).

At this point, I haven't really made any progress. The only tool I see for difference equations (diff_rec2) is designed for symbolic solutions of such systems, but in general economic models are not expressible in closed form. Runge-Kutta (rk) doesn't accept algebraic equations (as far as I see), and I'm not sure where to look next.

Ultimately, I would think this would be very straightforward given the direct forward-time computation nature of this and similar models. That said, I do want to avoid performing manual manipulations or turning this into a special case. I'm particularly interested in a general solution method for systems of such equations, as I plan to implement more complex models in the future, such as the McKinnon(1997) open economy.

---

Edit:

Thanks to Robert's (accepted) Answer, here's a fully-working copy-paste solution for the above forward-time difference equation example:

Y[t] := A[t]*K[t]^(1/3)*Ly[t]^(2/3);
K[t] := K[t - 1] + s*Y[t - 1] - d*K[t - 1];
A[t] := A[t - 1] + z*A[t - 1]*La[t - 1];
Ly[t] := L - La[t];
La[t] := l*L;

s : 0.15;
d : 0.07;
z : 0.02;
l : 0.05;
L : 100.0;

A[0] : 100.0;
K[0] : 100.0;

sol : makelist ([Y[n], K[n], A[n], Ly[n], La[t]], n, 0, 30);
v : makelist ([p-1, sol[p][1]], p, 1, 30);

plot2d ([discrete,v,0,30], logy);

Answer 1

Well, if the goal is to compute a solution for the equations given, I think you can do that with so-called memoizing functions (i.e., functions which compute a result and remember it). In Maxima such functions are defined by f[k] := ... (instead of f(k) := ... for ordinary, non-memoizing functions). In this case I think you would have:

Y[t] := A[t]*K[t]^(1/3)*Ly[t]^(2/3);
K[t] := K[t - 1] + s*Y[t - 1] - d*K[t - 1];
A[t] := A[t - 1] + z*A[t - 1]*La[t - 1];
Ly[t] := L - La[t];
La[t] := l*L;

 Y[0] : 9663.8253;
 K[0] : 100;
 A[0] : 100;
Ly[0] : 95.0;
La[0] : 5.0;

s: 0.15;
d: 0.07;
z: 0.02;
l: 0.05;
L: 100.0;

and then you can compute e.g. makelist ([Y[n], K[n], A[n], Ly[n], La[t]], n, 1, 10);

That works if you can rearrange the equations yourself; in this case it's easy. If you need to solve other problems and that approach fails, you'll need something more powerful. I don't know what to suggest. My advice is to ask on the maxima-discuss@lists.sourceforge.net mailing list where most discussion about Maxima takes place.