In[]:=
ClearAll["Global`*"];​​u[x_]:=(x^(1-gamma))/(1-gamma);​​du[x_]:=D[u[xx],xx]//.{xx->x};​​d=1;m=100;gamma=1.5;p=2;​​P1=Table[{p*i,p,p/i},{i,{1.2}}]//Flatten;​​B1=Join[{2*ConstantArray[1,9]},{KroneckerProduct[P1,ConstantArray[1,3]]//Flatten},{Table[{P1*i,P1,P1/i},{i,{1.2}}]//Flatten}];​​B1t=Transpose[B1];​​B2=Join[{a*ConstantArray[1,9]},{KroneckerProduct[Table[Symbol["b"<>ToString[i]],{i,1,3,1}],ConstantArray[1,3]]//Flatten},{Table[Symbol["c"<>ToString[i]],{i,1,9,1}]//Flatten}];​​B2t=Transpose[B2];​​BBig=Thread[Dot[B1t,B2t]];​​MBig=m*ConstantArray[1,9];​​cvars=B2[[3]];​​sol=Solve[BBigMBig,cvars];​​cvars=cvars/.sol//Flatten;​​cpath=Table[{Symbol["b"<>ToString[j]],Table[Symbol["c"<>ToString[i]],{i,3j-2,3j,1}]}//Flatten,{j,1,3,1}];​​Egen[y_,yH_,yM_,yL_]:=u[y]+(d/3)(u[yH]+u[yM]+u[yL]);​​Lpd=Egen@@Thread[cpath];​​Lpd=Lpd/.sol//Flatten;​​{Eqb1,Eqb2,Eqb3}=Lpd;​​DEqa=du[a]/P1[[2]]+(d/3)(du[b1]/P1[[1]]+du[b2]/P1[[2]]+du[b3]/P1[[3]])0;​​DEqb1=D[Eqb1,b1]0;​​DEqb2=D[Eqb2,b2]0;​​DEqb3=D[Eqb3,b3]0;​​
In[]:=
FindRoot[{DEqa,DEqb1,DEqb2,DEqb3},{a,5},{b1,5},{b2,5},{b3,5},MaxIterations->500]
FindRoot
:The search direction -5.22016×
107
10
+2.88475×
107
10
,2.17955×
107
10
-1.1571×
107
10
,2.61983×
107
10
-1.33261×
107
10
,3.14185×
107
10
-1.48373×
107
10
 is not a descent direction for the merit function. The step will be taken without the line search.
FindRoot
:The search direction -8.70027×
107
10
+4.80791×
107
10
,3.63259×
107
10
-1.92849×
107
10
,4.36638×
107
10
-2.22101×
107
10
,5.23642×
107
10
-2.47288×
107
10
 is not a descent direction for the merit function. The step will be taken without the line search.
FindRoot
:The search direction -1.45004×
108
10
+8.01319×
107
10
,6.05431×
107
10
-3.21415×
107
10
,7.27729×
107
10
-3.70169×
107
10
,8.72736×
107
10
-4.12146×
107
10
 is not a descent direction for the merit function. The step will be taken without the line search.
General
:Further output of FindRoot::bddir will be suppressed during this calculation.
FindRoot
:Failed to converge to the requested accuracy or precision within 500 iterations.
Out[]=
{a-1.66535×
111
10
+9.20302×
110
10
,b16.95328×
110
10
-3.69141×
110
10
,b28.35785×
110
10
-4.25133×
110
10
,b31.00232×
111
10
-4.73343×
110
10
}
In[]:=
BigU=u[a]+d/3(u[b1]+u[b2]+u[b3])+d^2/9(u[c1]+u[c2]+u[c3]+u[c4]+u[c5]+u[c6]+u[c7]+u[c8]+u[c9])/.sol;​​fun[ai_]:=(​​sol2=NSolve[{DEqb1/.aai,DEqb2/.aai,DEqb3/.aai},{b1,b2,b3}];​​Ufin=BigU/.({sol2,aai}//Flatten);​​Ufin//Flatten​​)​​FindMaximum[fun[as],{as,10}]
FindMaximum
:The function value 1.18246+0.000146503 is not a real number at {as} = {10.}.
Out[]=
FindMaximum[fun[as],{as,10}]