API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 26th, 2019
102
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 6.23 KB | None | 0 0
raw download clone embed print report
  1. jacob=D[system,{vars}]; (*for the Jacobian matrix*)
  2. Inverse[jacob]; (*inverse of the Jacobian*)
  3. Det[Jacob]; (*determinant of Jacobian*)
  4.  
  5. FindRoot::jsing: Encountered a singular Jacobian at the point` {EG,EXR,FSAV,IADJ,PA1,PA2,PD1,PD2,PE1,PM2,PQ1,PQ2,PVA1,PVA2,PX1,PX2,QA1,QA2,QD1,QD2,QE1,QF11,QF12,QF21,QF22,QH11,QH12,QH21,QH22,QINT11,QINT12,QINT21,QINT22,QINV1,QINV2,QM2,QQ1,QQ2,QX1,QX2,WALRAS,WF1,YF11,YF12,YF21,YF22,YG,YH1,YH2} = {109.,1.,4.,1.,1.,1.,1.,1.,1.,1.,1.,1.,0.6,0.5,1.,1.,279.,394.,249.,394.,30.,72.,105.,73.,135.,28.,47.,157.,87.,84.,55.,50.,99.,28.,85.,105.,259.,538.,279.,394.,0.,1.,95.,125.,82.,83.,109.,285.,186.}. Try perturbing the initial point(s).`
  6.  
  7. system={PM2 - EXR pwm2 (1 + tm2), PE1 - EXR pwe1 (1 - te1),
  8. PQ1 QQ1 - PD1 QD1 (1 + tq1),
  9. PQ2 QQ2 - (PD2 QD2 + PM2 QM2) (1 + tq2), -PD1 QD1 - PE1 QE1 +
  10. PX1 QX1, -PD2 QD2 + PX2 QX2, PA1 - PX1 [Theta]11 - PX2 [Theta]12,
  11. PA2 - PX1 [Theta]21 - PX2 [Theta]22, -PA1 + ica11 PQ1 + ica21 PQ2 +
  12. PVA1, -PA2 + ica12 PQ1 + ica22 PQ2 + PVA2,
  13. QA1 - ad1 QF11^[Alpha]11 QF21^[Alpha]21,
  14. QA2 - ad2 QF12^[Alpha]12 QF22^[Alpha]22,
  15. WF1 WFDIST11 - (PVA1 QA1 [Alpha]11)/QF11,
  16. WF1 WFDIST12 - (PVA2 QA2 [Alpha]12)/QF12,
  17. WF2 WFDIST21 - (PVA1 QA1 [Alpha]21)/QF21,
  18. WF2 WFDIST22 - (PVA2 QA2 [Alpha]22)/QF22, -ica11 QA1 +
  19. QINT11, -ica12 QA2 + QINT12, -ica21 QA1 + QINT21, -ica22 QA2 +
  20. QINT22, QX1 - QA1 [Theta]11 - QA2 [Theta]21,
  21. QX2 - QA1 [Theta]12 - QA2 [Theta]22,
  22. QQ2 - aq2 (QD2^-[Rho]q2 (1 - [Delta]q2) +
  23. QM2^-[Rho]q2 [Delta]q2)^(-1/[Rho]q2),
  24. QM2/QD2 - ((PD2 [Delta]q2)/(PM2 (1-[Delta]q2)))^(1/(1 + [Rho]q2)),
  25. -QD1 + QQ1,
  26. QX1 - at1 (QD1^[Rho]t1 (1 - [Delta]t1) + QE1^[Rho]t1 [Delta]t1)^(1/
  27. [Rho] t1),
  28. QE1/QD1 - ((PE1 (1 - [Delta]t1))/(PD1 [Delta]t1))^(
  29. 1/(-1 + [Rho]t1)), -QD2 +
  30. QX2, -shry11 (QF11 WF1 WFDIST11 + QF12 WF1 WFDIST12) +
  31. YF11, -shry12 (QF21 WF2 WFDIST21 + QF22 WF2 WFDIST22) +
  32. YF12, -shry21 (QF11 WF1 WFDIST11 + QF12 WF1 WFDIST12) +
  33. YF21, -shry22 (QF21 WF2 WFDIST21 + QF22 WF2 WFDIST22) +
  34. YF22, -tr13 - EXR tr14 - YF11 - YF12 + YH1, -tr23 - EXR tr24 -
  35. YF21 - YF22 + YH2, QH11 - ((1 - mps1) (1 - ty1) YH1 [Beta]11)/PQ1,
  36. QH12 - ((1 - mps2) (1 - ty2) YH2 [Beta]12)/PQ1,
  37. QH21 - ((1 - mps1) (1 - ty1) YH1 [Beta]21)/PQ2,
  38. QH22 - ((1 - mps2) (1 - ty2) YH2 [Beta]22)/PQ2, -IADJ qinv1 +
  39. QINV1, -IADJ qinv2 + QINV2, -EXR pwe1 QE1 te1 - EXR pwm2 QM2 tm2 -
  40. PD1 QD1 tq1 - (PD2 QD2 + PM2 QM2) tq2 - EXR tr34 + YG - ty1 YH1 -
  41. ty2 YH2, EG - PQ1 qg1 - PQ2 qg2 - tr13 - tr23, QF11 + QF12 - QFS1,
  42. QF21 + QF22 - QFS2, -qg1 - QH11 - QH12 - QINT11 - QINT12 - QINV1 +
  43. QQ1, -qg2 - QH21 - QH22 - QINT21 - QINT22 - QINV2 + QQ2,
  44. FSAV + pwe1 QE1 - pwm2 QM2 + tr14 + tr24 + tr34, -EG + EXR FSAV -
  45. PQ1 QINV1 - PQ2 QINV2 - WALRAS + YG + mps1 (1 - ty1) YH1 +
  46. mps2 (1 - ty2) YH2, -cpi + cwts1 PQ1 + cwts2 PQ2};
  47.  
  48. vars = {EG, EXR, FSAV, IADJ, PA1, PA2, PD1, PD2, PE1, PM2, PQ1, PQ2,
  49. PVA1, PVA2, PX1, PX2, QA1, QA2, QD1, QD2, QE1, QF11, QF12, QF21,
  50. QF22, QH11, QH12, QH21, QH22, QINT11, QINT12, QINT21, QINT22, QINV1,
  51. QINV2, QM2, QQ1, QQ2, QX1, QX2, WALRAS, WF1, YF11, YF12, YF21,
  52. YF22, YG, YH1, YH2};
  53.  
  54. parameters={[Rho]t1 -> 0.96, [Rho]t2 -> 1.1, [Rho]q1 -> 0.43,
  55. [Rho]q2 -> 0.43, [Alpha]11 -> 0.496552, [Alpha]21 ->
  56. 0.503448, [Alpha]12 -> 0.4375, [Alpha]22 -> 0.5625, [Beta]11 ->
  57. 0.153846, [Beta]21 -> 0.846154, [Beta]12 -> 0.347518, [Beta]22 ->
  58. 0.652482, [Theta]11 -> 1., [Theta]12 -> 0., [Theta]21 ->
  59. 0., [Theta]22 -> 1., [Delta]t1 -> 0.742333, [Delta]q2 ->
  60. 0.340471, shry11 -> 0.536723, shry21 -> 0.463277,
  61. shry12 -> 0.600962, shry22 -> 0.399038, ad1 -> 3.84818,
  62. ad2 -> 3.25772, at1 -> 2.56505, aq2 -> 2.42675, ica11 -> 0.301075,
  63. ica21 -> 0.179211, ica12 -> 0.139594, ica22 -> 0.251269,
  64. cwts1 -> 0.235119, cwts2 -> 0.764881, mps1 -> 0.264151,
  65. mps2 -> 0.220994, tq1 -> 0.0401606, tq2 -> 0.0371747,
  66. ty1 -> 0.0701754, ty2 -> 0.0268817, te1 -> 0, tm2 -> 0.371429,
  67. qg1 -> 13, qg2 -> 67, qinv1 -> 28, qinv2 -> 85, tr13 -> 25,
  68. tr23 -> 5, tr14 -> 40, tr24 -> 16, tr34 -> 15, QFS1 -> 177,
  69. QFS2 -> 208, WFDIST11 -> 0.711, WFDIST12 -> 1.579, WF2 -> 1,
  70. WFDIST21 -> 1, WFDIST22 -> 1, cpi -> 1.036, pwe1 -> 1, pwm2 -> 1};
  71.  
  72. givenPoint= {
  73. {EG, 109}, {EXR, 1}, {FSAV, 4}, {IADJ, 1}, {PA1, 1}, {PA2,
  74. 1}, {PD1, 1}, {PD2, 1}, {PE1, 1}, {PM2, 1}, {PQ1, 1}, {PQ2,
  75. 1}, {PVA1, 0.6}, {PVA2, 0.5}, {PX1, 1}, {PX2, 1}, {QA1,
  76. 279}, {QA2, 394}, {QD1, 249}, {QD2, 394}, {QE1, 30}, {QF11,
  77. 72}, {QF12, 105}, {QF21, 73}, {QF22, 135}, {QH11, 28}, {QH12,
  78. 47}, {QH21, 157}, {QH22, 87}, {QINT11, 84}, {QINT12, 55}, {QINT21,
  79. 50}, {QINT22, 99}, {QINV1, 28}, {QINV2, 85}, {QM2, 105}, {QQ1,
  80. 259}, {QQ2, 538}, {QX1, 279}, {QX2, 394}, {WALRAS, 0}, {WF1,
  81. 1}, {YF11, 95}, {YF12, 125}, {YF21, 82}, {YF22, 83}, {YG,
  82. 109}, {YH1, 285}, {YH2, 186}
  83. };
  84.  
  85. FindRoot[system/.parameters, givenPoint]
  86.  
  87. LinearSolve[
  88. jacob /. parameters /. Rule @@@ initialVal]["ConditionNumber"]
  89.  
  90. (* 4.75451*10^16 <-- condition number > 10.^MachinePrecision *)
  91.  
  92. sol = FindRoot[SetPrecision[system /. parameters, 24],
  93. SetPrecision[givenPoint, 24], PrecisionGoal -> 8,
  94. WorkingPrecision -> 24] // N
  95. (*
  96. {EG -> 116.757, EXR -> 1.57378, FSAV -> 10.052, IADJ -> 1.51676,
  97. PA1 -> 0.796624, PA2 -> 1.09585, PD1 -> 0.505239, PD2 -> 1.09585,
  98. PE1 -> 1.57378, PM2 -> 2.15832, PQ1 -> 0.52553, PQ2 -> 1.19291,
  99. PVA1 -> 0.424617, PVA2 -> 0.722743, PX1 -> 0.796624, PX2 -> 1.09585,
  100. QA1 -> 330.749, QA2 -> 337.713, QD1 -> 361.553, QD2 -> 337.713,
  101. QE1 -> 51.3489, QF11 -> 104.764, QF12 -> 72.2356, QF21 -> 70.705,
  102. QF22 -> 137.295, QH11 -> 61.6308, QH12 -> 97.7302, QH21 -> 149.331,
  103. QH22 -> 80.8366, QINT11 -> 99.5802, QINT12 -> 47.1428,
  104. QINT21 -> 59.2738, QINT22 -> 84.8569, QINV1 -> 42.4692,
  105. QINV2 -> 128.924, QM2 -> 132.401, QQ1 -> 361.553, QQ2 -> 570.222,
  106. QX1 -> 330.749, QX2 -> 337.713, WALRAS -> -5.72072*10^-14,
  107. WF1 -> 0.936218, YF11 -> 94.7432, YF12 -> 125., YF21 -> 81.7784,
  108. YF22 -> 82.9999, YG -> 159.551, YH1 -> 307.694, YH2 -> 194.959}
  109. *)
  110.  
  111. system /. parameters /. sol
  112. (*
  113. {0., 0., 0., 0., 0., 0., 0., 0., 0., 2.22045*10^-16, -5.68434*10^-14, 5.68434*10^-14,
  114. 1.11022*10^-16, 0., 0., -2.22045*10^-16, 0., 0., 0.,
  115. -1.42109*10^-14, 0., 0., 1.13687*10^-13, 1.11022*10^-16, 0., -5.68434*10^-14,
  116. 2.77556*10^-17, 0., -1.42109*10^-14, 0., -1.42109*10^-14, 0., 0.,
  117. 1.42109*10^-14, 1.42109*10^-14, 0., 0., 0., 0., 0., 0., 0., 0., -2.84217*10^-14, 0.,
  118. 5.68434*10^-14, 0., 0., 0.}
  119. *)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!