/* Econ2.mac Maxima software for Economic Analysis Ted Woollett, April 7, 2022 http://home.csulb.edu/~woollett/ http://home.csulb.edu/~woollett/eam.html Copyright (C) 2021, 2022, Edwin L. Woollett This program is free software: you can redistribute it and/or modify it under the terms of the GNU GENERAL PUBLIC LICENSE, Version 2, June 1991, as published by the Free Software Foundation. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.fsf.org/licensing/. */ /* DE1 (a,b,t,y0) or DE1 (a,b,t) returns the solution of the first order difference equation y[t] = b*y[t-1] + a. The optional fourth argument y0 is taken as a given y[0]. If y0 is not specified, the solution is returned in a form which includes an undefined constant %A. This function is used in Dowling17.wxmx. */ DE1 ( [argL] ) := block ([ %r ], if length (argL) = 4 then ( if argL[2] = 1 then argL[4] + argL[1]*argL[3] else ( %r : argL[1]/(1 - argL[2]), (argL[4] - %r)*argL[2]^argL[3] + %r)) else if length (argL) = 3 then ( kill (%A), if argL[2] = 1 then %A + argL[1]*argL[3] else ( %r : argL[1]/(1 - argL[2]), %A*argL[2]^argL[3] + %r)) else ( print ("syntax error"), done))\$ /* linear second order difference equations with constant coefficients */ /* Given the 2nd order linear difference equation y[t] + b1*y[t-1] + b2*y[t-2] = a, with constants b1,b2, a, yp[t] is returned by the Maxima function ypart (b1,b2,a,t). */ ypart (_b1,_b2,_a, _t) := (if _a = 0 then 0 else if (_b1 + _b2 = -1 and _b1 = -2) then (_a/2)*_t^2 else if (_b1 + _b2 = -1) then (_a*_t)/(2 + _b1) else _a/(1 + _b1 + _b2))\$ /* and yc[t] is given by the Maxima function ycompl (b1, b2, t). */ ycompl (_b1,_b2,_t) := block ( [%r1,%r2,%r,sarg,%g,%h,%k,%th], kill(%A1,%A2), sarg : _b1^2 - 4*_b2, /* assume coefficients are numerical here */ if sarg = 0 or float (sarg) = 0.0 then ( /* case repeated real roots */ %r : -_b1/2, %A1*%r^_t + %A2*_t*%r^_t) else if sarg > 0 then ( /* case distinct real roots */ %r1 : (-_b1 + sqrt(sarg))/2, %r2 : (-_b1 - sqrt(sarg))/2, %A1*%r1^_t + %A2*%r2^_t) else ( /* case sarg < 0, complex roots */ %g : -_b1/2, %h : sqrt(-sarg)/2, %k : sqrt (%g^2 + %h^2), if %g = 0 then %th : asin(%h/%k) else %th : acos (%g/%k), (%k^_t)*(%A1*cos (%th*_t) + %A2*sin (%th*_t))))\$ /* indefinite solution ysoln (b1,b2,a,t) given the 2nd order linear difference equation y[t] + b1*y[t-1] + b2*y[t-2] = a, */ ysoln (bb1,bb2,aa,tt) := (ypart (bb1,bb2,aa,tt) + ycompl (bb1,bb2,tt))\$ /* do loop code to produce a list of values y[t]: ytList(b1,b2,a,y0,y1,tmax), given the 2nd order linear difference equation y[t] + b1*y[t-1] + b2*y[t-2] = a */ ytlist(bb1,bb2,aa,yy0,yy1, ttmax) := block ([ ytm1 : yy1, ytm2 : yy0, yyt,yyp, LL ], if freeof(t, aa) then ( yyp : ypart (bb1,bb2,aa,t), print (" ye = ",yyp)), LL : [yy1, yy0], ytm1 : yy1, ytm2 : yy0, for tt:2 thru ttmax do( yyt : -bb1*ytm1 - bb2*ytm2 + subst (tt,t,aa), LL : cons (yyt, LL), ytm2 : ytm1, ytm1 : yyt), reverse (LL))\$ /* DE1S (A, B, t, Y0) or DE1S (A, B, t) returns Y[t], the matrix column vector solution of a set of n first order difference equations expressed in matrix form as Y[t] = A . Y[t-1] + B. Y0 is the optional matrix column vector of initial values. Y[t] has with n components if A is a square n x n matrix, and is returned in definite form if the fourth arg Y0 is included, else in terms of n constants _k[j]. */ DE1S ( [argL] ) := block ([ _A, _B, _t, _n, _Ye, _eval, _evec, _Yindef], local(_r,_V,_k), _A : argL[1], _B : argL[2], _t : argL[3], _n : length (_A), _Ye : invert (ident(_n) - _A) . _B, /* display (_Ye), */ [_eval, _evec] : eigenvectors (_A), if length (_eval[2]) < _n then return ("DE1S is not able to handle repeated roots"), /* display (_eval, _evec), */ for j thru _n do ( _r[j] : _eval[1][j], _V[j] : cvec (_evec[j][1])), _Yindef : _Ye + sum (_k[j]*_r[j]^_t*_V[j], j, 1, _n), /* display (_Yindef), */ if length (argL) = 3 then _Yindef else at (_Yindef, colVecSolve (at (_Yindef, _t = 0), argL[4])))\$ /* DE2S (A1, A2, B, t, Y0) The Maxima function DE2S (A1, A2, B, t, Y0), uses matrix methods to solve for the solution of A1 . Y[t] = A2 . Y[t - 1] + B, in which Y[t] is a matrix column vector with n elements depending on t, A1 and A2 are square n x n matrices of numerical elements, and B is either a given n element matrix column vector (with constant numerical values), or the number 0. Finally Y0 is an n element numerical matrix column vector giving the desired initial values Y(0). For a problem in which there is no B term, and we are solving A1 . Y[t] = A2 . Y[t-1], with Y[0] = Y0, we can either use DE2S (A1,A2, 0, Y0) or DE2S (A1,A2, zeromatrix (n, 1), Y0), in which n = length(A) = length (Y0). The solution method is to multiply every term by invert(A1) (provided det(A1) # 0), which reduces the set of equations to the form Y[t] = D . Y[t-1] + E, and then call DE1S. */ DE2S (_A1, _A2, _B, _t, _Y0) := block ( if det (_A1) = 0 then return (" DE2S method depends on det (A1) not being zero"), DE1S (invert (_A1) . _A2, invert (_A1) . _B, _t, _Y0))\$ /* Ytlist (A, B, Y0, tmax) produces a column vector whose jth component is the list of values [yj[0], yj[1], yj[2],...,yj[tmax] ] using directly the first order difference equation Y[t] = A . Y[t-1] + B repeatedly, starting with Y[1] = A . Y0 + B, without using the eigenvalues or eigenvectors of A. */ Ytlist (AA, BB, YY0, ttmax) := block ([nn,YYp,YYt, CCL : [] ], local (yyL), nn : length (AA), for j thru nn do yyL[j] : [YY0[j,1]], YYp : YY0, for tv thru ttmax do ( YYt : AA . YYp + BB, for j thru nn do yyL[j] : cons (YYt[j,1], yyL[j] ), YYp : YYt), for j thru nn do yyL[j] : [reverse ( yyL[j] )], for j thru nn do CCL : cons (yyL[j], CCL), apply ('matrix, reverse (CCL) ) )\$ /* pIUlist (g, j, βk, m, y0, y1, tmax), produces a list of values of the inflation rate p[t], p[0] = y0, p[1] = y1, p[2], ..., p[tmax], corresponding to given numerical values of g, j, β*k, m, y0, y1, and tmax. This is based on the model of inflation and unemployment in discrete time (Chiang and Wainwright, Sec 18.3) discussed in the worksheet Dowling18C.wxmx. After calculating b1, b2, and a, corresponding to the 2nd order linear difference equation p[t] + b1*p[t-1] + b2*p[t-2] = a or equivalently p[t+2] + b1*p[t+1] + b2*p[t] = a, this function calls ytlist (b1,b2,a,y0,y1,tmax). */ pIUlist (gg,jj,betak,mm,%y0,%y1,%tmax) := block ([bb1,bb2,aa], /* check params */ if betak <= 0 then return("need β*k > 0"), if gg < 0 or gg > 1 then return("need 0 < g <= 1"), if jj < 0 or jj > 1 then return("need 0 < j <= 1"), if not numberp (mm) then return (" m must be a number"), aa : jj*mm*betak/(1 + betak), bb1 : -(1 + gg*jj + (1 - jj)*(1 + betak))/(1 + betak), bb2 : (1 - jj*(1 - gg))/(1 + betak), ytlist (bb1,bb2,aa, %y0, %y1,%tmax))\$ /* Multiplier-Accelerator model do loop code MAlist (, b, IaG0, Y0, Y1, tmax) produces a list of lists of the form [t, Ct, Idt, IaG0, Yt] */ MAlist (aa,bb,iag0,yy0,yy1, ttmax) := block ([ ytm1 : yy1, ytm2 : yy0, Ct, Idt, Yt, LL ], Ye : iag0/(1 - bb), display (Ye), LL : [ ["t", "Ct","Idt","Ia + G0","Yt"] ], LL : cons ([0, 0, 0, iag0, yy0], LL), LL : cons ([1, bb*yy0, 0, iag0, yy1], LL), ytm1 : yy1, ytm2 : yy0, for tt:2 thru ttmax do( Ct : bb*ytm1, Idt : aa*bb*(ytm1 - ytm2), Yt : Ct + Idt + iag0, LL : cons ([tt, Ct, Idt,iag0, Yt], LL), ytm2 : ytm1, ytm1 : Yt), reverse (LL))\$ /* if maL is the list produced by MAlist, apply ('matrix, maL) will turn the output of MAlist into a formatted table. */ /* makeLists(maL) produces sublists tL, CtL, IdtL, IaG0L, YtL from the output of MAlist. */ aL(LL, num) := makelist (LL[j][num], j, 2, length (LL))\$ makeLists (alist) := (tL : aL(alist,1), CtL : aL(alist,2), IdtL : aL(alist,3), IaG0L : aL(alist,4), YtL : aL(alist,5), done)\$ /* For 2nd order ODE's with constant coefficients and terms with the general form y'' + A*y' + B*y(t) = C,an alternative approach to using ode2 is to design a Maxima function which requires us to decide (ahead of time) on the relative sizes of A,B,C (or express them symbolically) and ask for the specific cases: complex roots (sines and cosines) if A^2 < 4*B, distinct real roots (if A^2 > 4*B), or double real roots (if A^2 = 4*B). For the arg 'type' in Lode2 then use either complex, real, or double. This function is used in Dowling18A.wxmx. */ /* linear ode2 with constant coefficients y'' + A y' + B y(t) = C Lode2(y,t,type,A,B,C) where type = either real, complex, or double. */ Lode2(%y,%t,%type,%A,%B,%C) := block ([_omega, _r,_term, _arg ], if %A = 0 or %B = 0 then ode2('diff(%y,%t,2) + %A*'diff(%y,%t) + %B = %C, %y, %t) else ( kill(%k1,%k2), _r : %A/2, if %C = 0 then _term : 0 else ( _term : %C/%B, _term : (-num(_term))/(-denom(_term))), if %type = real then ( print(" this assumes that ", %A^2, " > ",4*%B ), _arg : %A^2 - 4*%B, if numberp (_arg) then if _arg < 0 then return ("assumption incorrect"), /* print ("_arg = ", _arg), */ if dofactor then _arg : factorsum (_arg), /* print ("_arg = ", _arg), */ _omega : sqrt(_arg)/2, /* print ("_omega = ", _omega), */ if doratsimp then %y = %k1*exp(- ratsimp((_r - _omega))*%t) + %k2*exp(- ratsimp((_r + _omega))*%t) + _term else %y = %k1*exp(- (_r - _omega)*%t) + %k2*exp(- (_r + _omega)*%t) + _term) else if %type = complex then ( print (" this assumes that ", 4*%B, " > ", %A^2), _arg : 4*%B - %A^2, if numberp (_arg) then if _arg < 0 then return ("assumption incorrect"), /* print ("_arg = ", _arg), */ if dofactor then _arg : factorsum (_arg), /* print ("_arg = ", _arg), */ _omega : sqrt (_arg)/2, /* print ("_omega = ", _omega), */ %y = exp(-_r*%t)*(%k1*sin(_omega*%t) + %k2*cos(_omega*%t)) + _term) else if %type = double then ( print (" this assumes that ", %A^2, " = ", 4*%B), if %C = 0 then _term : 0 else ( _term : 4*%C/%A^2, _term : (-num(_term))/(-denom(_term))), %y = (%k1 + %k2*%t)*exp (- _r*%t) + _term ) else print ("type should be real, complex, or double")))\$ /* croots (A,B) 1. returns the solutions of the equation r^2 + A*r + B = 0. 2. returns the characteristic roots of the linear second order ordinary differential equation (ode) with constant coefficients: y'' + A*y' + B*y = C. */ croots(%A,%B) := reverse (map ('rhs, solve (rr^2 + %A*rr + %B = 0, rr)))\$ /* with symbolic A,B: (%i6) croots (A,B); (%o6) [ (sqrt(A^2 - 4*B) -A)/2, -(sqrt(A^2 - 4*B)+A)/2 ] croots (A,B), numer; can also be used with numerical values of A and B. [r1, r2] : croots(A,B) assigns the first element of the list returned to the symbol r1, the second to r2. */ squarep (MM) := block (if length (MM) = length (transpose (MM)) then true else false)\$ det (MM) := block(if squarep (MM) then determinant (MM) else return ("matrix must be square"))\$ mtrace (MM) := block ( if squarep (MM) then sum (MM[j, j], j,1, length (MM)) else return (" matrix must be a square matrix"))\$ /* Mcroots(M) assumes length(M) = 2 */ Mcroots (MM) := block ([MSqrt, MTr], if squarep (MM) then ( MTr : mtrace (MM), MSqrt : sqrt ( MTr^2 - 4* abs (det (MM))), [ (MTr + MSqrt)/2, (MTr - MSqrt)/2]) else return (" matrix must be a square matrix"))\$ /* following functions also defined in mbe5.mac */ cvec(zzL) := ( if not listp(zzL) then ( print(" arg of cvec should be a list"), return()), transpose (matrix (zzL) ))\$ alias(lme, list_matrix_entries)\$ /* to_solve(Ab,xL) generate list of n equations from augmented matrix Ab with n rows and n+1 columns, given n-element list xL (%i26) xL; (%o26) [x1,x2] (%i27) Cd; (%o27) matrix([1,1,2],[0,1,1]) (%i28) to_solve(Cd, xL); (%o28) [x2+x1 = 2,x2 = 1] */ to_solve(Ab1, xL1) := block( [rL, eqn1, eqn_list : [], j], for j thru length(xL1) do ( rL : lme (row (Ab1,j)), eqn1 : rest(rL,-1) . xL1 = last(rL), eqn_list : cons (eqn1, eqn_list)), reverse (eqn_list))\$ /* solve_aug (M,L) solves for the n unknowns in a list L, given the augmented matrix M for a system of n equations. Calls to_solve defined above and then solve. */ solve_aug(Bc1,zL1) := block([ssL], ssL : to_solve(Bc1,zL1), ssL : solve(ssL, zL1), if length(ssL) = 1 then first(ssL) else ssL)\$ /* (%i31) xL; (%o31) [x1,x2] (%i32) Cd; (%o32) matrix([1,1,2],[0,1,1]) (%i33) solve_aug(Cd,xL); (%o33) [[x1 = 1,x2 = 1]] */ /* Solve for n unknown coefficients in n given equations implied by the matrix equation B = C, in which B and C are each column vectors of length n. */ colVecSolve(%B, %C) := block ([%eqnL : [], num : length (%B) ], if length (%C) # num then return (" column vecs must be same length"), for j thru num do %eqnL : cons (%B[j,1] = %C[j,1], %eqnL), solve (%eqnL)[1])\$ ODE1S_debug (%A, %B, %Y0) := block ([%ev,%mul,%evec,%n, %abort, %Yindef], local (%r, %V, %k), if not squarep (%A) then return (" A must be a square matrix"), %n : length (%A), if length (%Y0) # %n then return ("Y0 must be a matrix column vector of length ", %n), if numberp(%B) then if %B = 0 then %B : zeromatrix (%n, 1) else return ("B must be a zeromatrix of length ",%n," or the number 0"), if length (%B) # %n then return ("B must be a matrix column vector of length ",%n, " or else the number 0"), [%ev, %evec] : eigenvectors (%A), display (%ev, %evec), %mul : %ev[2], display (%mul), for j thru %n do ( print ("j = ", j," %mul[j] = ", %mul[j]), if %mul[j] # 1 then %abort : true), if %abort then return (" Found repeated eigenvalue"), for j thru %n do ( %r[j] : %ev[1][j], print ("j = ", j," r[j] = ",%r[j])), for j thru %n do ( %V[j] : cvec (%evec[j][1]), print ("j = ", j, " evec[j] = ", %V[j])), %Yindef : - invert (%A) . %B + sum (%k[j]*exp (%r[j]*t)*%V[j], j, 1, %n), display (%Yindef), %solns : colVecSolve (at (%Yindef, t = 0), %Y0), return ( at (%Yindef, %solns)))\$ /* ODE1S (A, B, Y0) The Maxima function ODE1S (A, B, Y0), uses matrix methods to solve for the solution of dY/dt = A . Y(t) + B, in which Y is a matrix column vector with n elements depending on t, A is a square n x n matrix of numerical elements, and B is either a given n element matrix column vector (with constant numerical values), or the number 0. Finally Y0 is an n element numerical matrix column vector giving the desired initial values Y(0). For a problem in which there is no B term, and we are solving dY/dt = A . Y(t), with Y(0) = Y0, we can either use ODE1S (A, 0, Y0) or ODE1S (A, zeromatrix (n, 1), Y0), in which n = length(A) = length (Y0). */ ODE1S (%A, %B, %Y0) := block ([%ev,%mul,%evec,%n, %abort : false, %Yindef], local (%r, %V, %k), if not squarep (%A) then return (" A must be a square matrix"), %n : length (%A), if length (%Y0) # %n then return ("Y0 must be a matrix column vector of length ", %n), if (matrixp(%B) and length (%B) # %n) then return ("B must be a matrix column vector of length ",%n," or the number 0"), if (numberp(%B) and %B = 0) or (matrixp(%B) and %B = zeromatrix (%n, 1)) then return (expand (matrixexp (%A, t) . %Y0)), /* case %B present in problem */ [%ev, %evec] : eigenvectors (%A), %mul : %ev[2], for j thru %n do if %mul[j] # 1 then %abort : true, if %abort then return (" Found repeated eigenvalue"), for j thru %n do ( %r[j] : %ev[1][j], %V[j] : cvec (%evec[j][1])), %Yindef : - invert (%A) . %B + sum (%k[j]*exp (%r[j]*t)*%V[j], j, 1, %n), %solns : colVecSolve (at (%Yindef, t = 0), %Y0), return ( at (%Yindef, %solns)))\$ /* (%i87) Ys : ODE1S (matrix ([5, -0.5], [-2, 5]), cvec ([-12, -24]), cvec ([12, 4])); (Ys) matrix([5*%e^(6*t)+4*%e^(4*t) + 3.0],[-10*%e^(6*t)+ 8*%e^(4*t) + 6.0]) (%i88) at (Ys, t = 0); (%o88) matrix([12.0],[4.0]) (%i89) Ys : ODE1S (matrix ([6,5], [1,2]), 0, cvec ([4, 1])); (Ys) matrix( [(25*%e^(7*t))/6 - %e^t/6],[(5*%e^(7*t))/6 + %e^t/6]) (%i90) at (Ys, t = 0); (%o90) matrix([4],[1]) (%i91) Ys : ODE1S (matrix ([6,5], [1,2]), zeromatrix (2,1), cvec ([4, 1])); (Ys) matrix([(25*%e^(7*t))/6 - %e^t/6],[(5*%e^(7*t))/6 + %e^t/6]) (%i92) at (Ys, t = 0); (%o92) matrix([4],[1]) */ /* ODE2S (A1, A2, B, Y0) The Maxima function ODE2S (A1, A2, B, Y0), uses matrix methods to solve for the solution of A1 . dY/dt = A2 . Y(t) + B, in which Y is a matrix column vector with n elements depending on t, A1 and A2 are square nxn matrices of numerical elements, and B is either a given n element matrix column vector (with constant numerical values), or the number 0. Finally Y0 is an n element numerical matrix column vector giving the desired initial values Y(0). For a problem in which there is no B term, and we are solving A1 . dY/dt = A2 . Y(t), with Y(0) = Y0, we can either use ODE2S (A1,A2, 0, Y0) or ODE2S (A1,A2, zeromatrix (n, 1), Y0), in which n = length(A) = length (Y0). The solution method is to multiply every term by invert(A1) (provided det(A1) # 0), which reduces the set of equations to the form dY/dt = D . Y + E, and then call ODE1S. */ ODE2S (_A1, _A2, _B, _Y0) := block ( if det (_A1) = 0 then return (" ODE2S method depends on det (A1) not being zero"), ODE1S (invert (_A1) . _A2, invert (_A1) . _B, _Y0))\$ /* fll(L) returns the first element of the list L, the last element of the list L, and the length of the list L */ fll(x) := [first(x),last(x),length(x)]\$ /* Dowling ch 20 */ NumSufficient(d11,d1,d21,d2) := block ( if not numberp (d11) then display(d11), if not numberp (d1) then display (d1), if not numberp (d21) then display(d21), if not numberp (d2) then display (d2), if (d11 < 0 and d1 > 0) and (d21 < 0 and d2 > 0) then print ("global maximum") else if (d11 <= 0 and d1 >= 0) and (d21 <= 0 and d2 >= 0) then print ("relative maximum") else if (d11 > 0 and d1 > 0) and (d21 > 0 and d2 > 0) then print ("global minimum") else if (d11 >= 0 and d1 >= 0) and (d21 >= 0 and d2 >= 0) then print ("relative minimum") else print ("neither maximized nor minimized: saddle point"))\$ Extremal(FF) := block ([ffx,ffxp,ffxpt,ffxpx,ffxpxp,xsoln,varOde, AA,BB, CCL : [ ], CC, Eqn], kill(%k1), ffx : diff (FF,x), ffxp : diff (FF, xp), ffxpt : diff (ffxp, t), ffxpx : diff (ffxp,x), if not lfreeof ([x, xp, xpp], listofvars (ffxpx)) then (print(" nonlinear ode"), return(done )), ffxpxp : ratsimp (diff (ffxp, xp)), if not lfreeof ([x, xp, xpp], listofvars (ffxpxp)) then (print(" nonlinear ode"), return(done )), if details then display (ffx, ffxp, ffxpt, ffxpx, ffxpxp), /* print ("Euler's Equation"), print (ffx - ffxpt - ffxpx*xp - ffxpxp*xpp," = 0"), */ /* no derivatives case */ if ffxpx = 0 and ffxpxp = 0 then ( xsoln : solve(ffx - ffxpt,x)[1], /* display (xsoln), */ return (xsoln)) /* no xpp term case : fxpxp = 0, fxpx # 0 */ else if ffxpxp = 0 then ( xsoln : x = %k1 + integrate ( (ffx - ffxpt)/ffxpx, t), /* display (xsoln), */ return (xsoln)), /* case we have an xpp term */ varOde : ffx - ffxpt - ffxpx*xp - ffxpxp*xpp, if details then display (varOde), varOde : expand (varOde/coeff (varOde,xpp)), if details then display (varOde), AA : coeff(varOde, xp), BB : coeff (varOde,x), if details then display (AA, BB), for j thru length (varOde) do if lfreeof ([x, xp, xpp], part (varOde, j)) then CCL : cons (- part(varOde,j), CCL), if details then display (CCL), CC : apply ("+", CCL), if details then display (CC), Eqn : xpp + AA*xp + BB*x = CC, print ("ode: ", Eqn), Eqn : 'diff (x,t,2) + AA*'diff (x,t) + BB*x = CC, ode2(Eqn, x, t))\$ NumSuffCond(F) := block ([H1, H2, D11,D1, D21, D2 ], H1 : hessian (F, [x,xp]), if details then display (H1), D11 : H1[1,1], if details then display (D11), D1 : float (determinant (H1)), if details then display (D1), H2 : hessian (F, [xp, x]), if details then display (H2), D21 : H2[1,1], if details then display (D21), D2 : float (determinant (H2)), if details then display (D2), NumSufficient (D11, D1, D21, D2), done)\$ NumDynamic (F) := block ([cdex], cdex : Extremal (F), print (" candidate extremal: ", cdex), NumSuffCond (F), done)\$ Dynamic (F) := block ([H1, H2, D11,D1, D21, D2, Fx, Fxp, Fxpt, Fxpx, Fxpxp, xp, xpp, ode, A, B, C, eqn, soln ], H1 : hessian (F, [x,xp]), /* display (H1), */ D11 : H1[1,1], /* display (D11), */ D1 : float (determinant (H1)), /* display (D1), */ H2 : hessian (F, [xp, x]), /* display (H2), */ D21 : H2[1,1], /* display (D21), */ D2 : float (determinant (H2)), /* display (D2), */ sufficient (D11, D1, D21, D2), Fx : diff (F,x), Fxp : diff (F, xp), Fxpt : diff (Fxp, t), Fxpx : diff (Fxp,x), Fxpxp : diff (Fxp, xp), /* display (Fx, Fxp, Fxpt, Fxpx, Fxpxp), */ ode : Fx - Fxpt - Fxpx*xp - Fxpxp*xpp, ode : expand (ode/coeff (ode,xpp)), A : coeff(ode, xp), B : coeff (ode,x), /* display (A, B), */ for j thru length (ode) do if numberp (part (ode, j)) then C : - part(ode, j), if not numberp(C) then C : 0, /* display (C), */ eqn : 'diff (x,t,2) + A*'diff (x,t) + B*x = C, /* display (eqn), */ soln : ode2(eqn, x, t), /* display (soln), */ print (" candidate extremal is "), print (" x(t) = ", rhs(soln)), done)\$ /* Dowling Ch. 21 code */ ConcaveTest (FF, xx, yy) := block ([ d1,d11, dd1, d2, d21, dd2, %simply : false], d1 : hessian (FF, [xx, yy]), display (d1), d11 : d1[1,1], display (d11), dd1 : determinant (d1), display (dd1), if d11 < 0 and dd1 > 0 then ( print("strictly concave"), return (done)), if d11 <= 0 and dd1 >= 0 then ( /* print("check simply"), */ d2 : hessian (FF, [yy, xx]), display (d2), d21 : d2[1,1], display (d21), dd2 : determinant (d2), display (dd2), if d21 <= 0 and dd2 >= 0 then %simply : true), if %simply then print ("simply concave") else print (" no sufficient conditions"), done)\$ Zindef (%A, %B ) := block ([%ev,%mul,%evec, %n, %abort : false ], local (%r, %V, %k), %n : length (%A), [%ev, %evec] : eigenvectors (%A), %mul : %ev[2], for j thru %n do if %mul[j] # 1 then %abort : true, if %abort then return (" Found repeated eigenvalue"), for j thru %n do ( %r[j] : %ev[1][j], %V[j] : cvec (%evec[j][1])), - invert (%A) . %B + sum (%k[j]*exp (%r[j]*t)*%V[j], j, 1, %n))\$ Zindef_float (%A, %B ) := block ([%ev,%mul,%evec, %n, %abort : false ], local (%r, %V, %k), /* display (%A, %B), */ %n : length (%A), [%ev, %evec] : eigenvectors (%A), %mul : %ev[2], for j thru %n do if %mul[j] # 1 then %abort : true, if %abort then return (" Found repeated eigenvalue"), %ev : float(%ev), %evec : float(%evec), for j thru %n do ( %r[j] : %ev[1][j], %V[j] : cvec (%evec[j][1])), - invert (%A) . %B + sum (%k[j]*exp (%r[j]*t)*%V[j], j, 1, %n))\$ dofactor : true\$ doratsimp : true\$ details : false\$ fpprintprec:5\$ ratprint:false\$