Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Solving this in mathematica?

• To: mathgroup at smc.vnet.net
• Subject: [mg128066] Re: Solving this in mathematica?
• From: Bob Hanlon <hanlonr357 at gmail.com>
• Date: Thu, 13 Sep 2012 03:38:04 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• Delivered-to: l-mathgroup@wolfram.com
• Delivered-to: mathgroup-newout@smc.vnet.net
• Delivered-to: mathgroup-newsend@smc.vnet.net
• References: <20120912065824.420E96824@smc.vnet.net>

rbar = 0.006236 // Rationalize;

(* rt=r_bar You cannot use _ in a name.
I assume the RHS should be rbar *)

rt = rbar; (* rt is not used *)
k =
 0.95 // Rationalize;
sigmar = 0.002 // Rationalize;
betazr = -0.00014 // Rationalize;
sigmaz = 0.4 // Rationalize;
pi = 0.99 // Rationalize; (* pi is not used *)
chi =
 0.05 // Rationalize; (* chi is not used *)
Cbar = -3.7 // Rationalize;

(* Use RSolve to solve your recurrence equations *)

Clear[alpha1, alpha2];
sol1 = {alpha1[n], alpha2[n]} /. RSolve[{
      alpha1[n] == alpha1[n - 1] + alpha2[n - 1],
      alpha2[n] == k (alpha2[n - 1]),
      alpha1[0] == 0, alpha2[0] == 1},
     {alpha1[n], alpha2[n]}, n][[1]];

alpha1[n_] = sol1[[1]];
alpha2[n_] = sol1[[2]];

Clear[sigma1sq, sigma12, sigma2sq];
sol2 = {sigma1sq[n], sigma12[n], sigma2sq[n]} /. RSolve[{
       sigma1sq[n] == sigma2sq[n - 1] + 2 sigma12[n - 1] + sigmaz^2,
       sigma12[n] == k (sigma12[n - 1]) + k (sigma2sq[n - 1]) + betazr,
       sigma2sq[n] == (k^2) (sigma2sq[n - 1]) + sigmar^2,
       sigma1sq[0] == 0, sigma12[0] == 0, sigma2sq[0] == 0},
      {sigma1sq[n], sigma12[n], sigma2sq[n]}, n][[1]] //
   FullSimplify;

sigma1sq[n_] = sol2[[1]];
sigma12[n_] = sol2[[2]];
sigma2sq[n_] = sol2[[3]];

(* You did not specify a starting value for psi[n], I have assumed that \
psi[0] = 0 *)

Clear[phi1, phi2, psi];
sol3 = {phi1[n], phi2[n], psi[n]} /. RSolve[{
       phi1[n] == phi1[n - 1] + phi2[n - 1] + (sigmaz^2)/2,
       phi2[n] == k (phi2[n - 1]) + (1 - k) (rbar),
       psi[n] == phi1[n] - sigma1sq[n]/2,
       phi1[0] == 0, phi2[0] == 0, psi[0] == 0},
      {phi1[n], phi2[n], psi[n]}, n][[1]] //
   FullSimplify;

phi1[n_] = sol3[[1]];
phi2[n_] = sol3[[2]];
psi[n_] = sol3[[3]];

B[h_, r_] = Exp[-alpha1[h]* r - psi[h]] //
    PiecewiseExpand // Simplify;

NMinimize[Abs[Exp[Cbar - beta] Sum[(Pi^x) B[x, r],
     {x, 1, 1000}] - 1], {r, beta},
 WorkingPrecision -> 35]

{1.0000000000000000000000000000000000, {r ->
   311.81988579682125561865976837161440,
  beta -> 158.37558965184462887837320791605232}}

This apears to indicate that there is no solution to your equation
since the minimum is not zero. Although this may be just a local
minimum and the function may need further exploration.


Bob Hanlon


On Wed, Sep 12, 2012 at 2:58 AM, Leon <leon.he88 at gmail.com> wrote:
> I have following code in mathematica:
>
> ------------------------------------------------------------
> rbar = 0.006236
> rt = r_bar
> k = 0.95
> sigmar = 0.002
> betazr = -0.00014
> sigmaz = 0.4
> pi = 0.99
> chi = 0.05
> Cbar = -3.7
>
> alpha1[n_] := alpha1[n] = alpha1[n - 1] + alpha2[n - 1]
> alpha2[n_] := alpha2[n] = k (alpha2[n - 1])
> sigma1sq[n_] :=
>  sigma1sq[n] = sigma2sq[n - 1] + 2 sigma12[n - 1] + sigmaz^2
> sigma12[n_] :=
>  sigma12[n] = k (sigma12[n - 1]) + k (sigma2sq[n - 1]) + betazr
> sigma2sq[n_] := sigma2sq[n] = (k^2) (sigma2sq[n - 1]) + sigmar^2
> phi1[n_] := phi1[n] = phi1[n - 1] + phi2[n - 1] + (0.5) (sigmaz^2)
> phi2[n_] := phi2[n] = k (phi2[n - 1]) + (1 - k) (rbar)
> psi[n_] := psi[n] = phi1[n] - (0.5) (sigma1sq[n])
>
> alpha1[0] = 0
> alpha2[0] = 1
> sigma1sq[0] = 0
> sigma12[0] = 0
> sigma2sq[0] = 0
> phi1[0] = 0
> phi2[0] = 0
>
> B[h_, r_] := Exp[(-alpha1[h]) (r) - psi[h]]
> Exp[Cbar - beta] Sum[(Pi^x) B[x, r], {x, 1, 1000}]
>
> ----------------------------------------------------------------------
>
> and I am wondering if it is possible to solve the last line such that I have "r" as a function of "beta", satisfying
>
> Exp[Cbar - beta] Sum[(Pi^x) B[x, r], {x, 1, 1000}] = 1
>
> Because ultimately, I would need to integrate a function J[r] over "beta", so if I don't have "r" as a function of "beta", I don't know how to do the integration of J[r].
>
> Thank you!! L.
>

• References:
  • Solving this in mathematica?
    • From: Leon <leon.he88@gmail.com>