3
$\begingroup$

I have two functions $P_1,Q_1$ which are both depending on the variables $r_{T_m}, r_{i_f}$, which are theirselves functions of $T_m, mi_f$. I am interested in computing the Jacobian

$$\frac{\partial (P_1,Q_1)}{\partial(r_{T_m}^2,r_{i_f}^2)}\Bigg|_{r_{Tm}=f(T_m) \\ r_{i_f}=g(mi_f)}.$$

To do that, I first define the function

powerD[f_, x_^(k_.)] := powerD[f, {x^k, 1}];
powerD[f_, {x_^(k_.), 0}] := f;
powerD[f_, vars__] := Fold[powerD, f, {vars}];
powerD[f_, {x_^(k_.), n_Integer?Positive}] := Det[Append[Table[(j!/i!) Binomial[k i, j] x^(k i - j), {i, n - 1}, {j, n}], Table[D[f, {x, j}], {j, n}]]]/(k x^(k - 1))^Binomial[n + 1, 2];

(because I want to compute the derivative wrt the power of $r_{T_m}, r_{mi_f}$), then I want to substitute $r_{T_m}=f(T_m)$ and $r_{i_f}=g(mi_f)$. However, something in the code fails at this point. I think it is because instead of first computing the partial derivatives $\frac{\partial P_1}{\partial r^2_{T_m}}$ etc. and then substituting $r_{T_m}=f(T_m)$, it tries to compute $\frac{\partial P_1}{\partial f^2(T_m)}$.

I attach here the code ($P_1,Q_1$ are actually functions of other parameters too, but I substitute these other parameters with the needed numerical values before computing the Jacobian).

The function $P_1(r_{T_m},mi_f)$ is given by

P1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif]=-1.12469*10^-8 (1.71055*10^13 - 112.657 rif^2 + 112.657 rTm^2 - 9.60432*10^-10 (-5.66361*10^20 - 3.52759*10^9 rif^2 + 3.52759*10^9 rTm^2 + Sqrt[(5.66361*10^20 + 3.52759*10^9 rif^2 - 3.52759*10^9 rTm^2)^2 - 3.26026*10^16 (3.2715*10^26 - 3.83021*10^15 rif^2 + 12691.7 rif^4 - 4.07533*10^15 rTm^2 - 25383.3 rif^2 rTm^2 + 12691.7 rTm^4)]))

and the function $Q_1(r_{T_m},mi_f)$ is given by

Q1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif]=6.13447*10^-17 (-5.66361*10^20 - 3.52759*10^9 rif^2 + 3.52759*10^9 rTm^2 + Sqrt[(5.66361*10^20 + 3.52759*10^9 rif^2 - 3.52759*10^9 rTm^2)^2 - 3.26026*10^16 (3.2715*10^26 - 3.83021*10^15 rif^2 + 12691.7 rif^4 - 4.07533*10^15 rTm^2 - 25383.3 rif^2 rTm^2 + 12691.7 rTm^4)])

I compute the elements of the (2x2) Jacobian matrix as

J11[rTm_, rif_] := FullSimplify[powerD[P1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif], rTm^2]]
J12[rTm_, rif_] := FullSimplify[powerD[P1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif], rif^2]]
J21[rTm_, rif_] := FullSimplify[powerD[Q1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif], rTm^2]]
J22[rTm_, rif_] := FullSimplify[powerD[Q1[230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif], rif^2]]

and then I build the full Jacobian matrix as

J[rTm_, rif_] := {{J11[rTm, rif], J12[rTm, rif]}, {J21[rTm, rif],J22[rTm, rif]}}

Up to here everything is okay, but then when I substitute $$r_{T_m}=f(T_m)=\sqrt\frac{(230\cdot1.73)^4 + 4\cdot(230\cdot1.73)^2\cdot0.2\cdot100\cdot3.14\cdot Tm}{4\cdot(0.2)^2}$$

and

$$r_{i_f}=g(mi_f)=\sqrt\frac{230\cdot1.73\cdot mif\cdot100\cdot3.14 }{(0.2)^2 + (100\cdot3.14\cdot 0.0072)^2}$$

with the command

J[Sqrt[((230 1.73)^4 + 4 (230 1.73)^2 0.2 100 3.14 Tm)/(4 (0.2)^2)], Sqrt[((230 1.73) mif 100 3.14)/((0.2)^2 + (100 3.14 0.0072)^2)]]

I get the errors "6.25 (2.5066610^10+3.9771110^7 Tm) is not a valid variable", "24254.6mif is not a valid variable."

I have just started today using the software, so I am sure it is a pretty easy mistake to be solved, maybe I should assign the values differently to the various Jacobian elements $J_{i,j}$.

$\endgroup$

1 Answer 1

1
$\begingroup$

Consider Q[r1, r2] and P[r1, r2] (I am using simplified notation). The Jacobian is defined by:

D[{Q[r1, r2], P[r1, r2]}, {{r1, r2}}]

However, you want derivatives wrt. r1^2 and r2^2. Toward this aim we introduce new variables :

w1 = r1^2;
w2 = r2^2;

and Q and P are written: Q[Sqrt[w1], Sqrt[w2]] and P[Sqrt[w1], Sqrt[w2]].

With this we now get the Jacobian:

D[{Q[Sqrt[w1], Sqrt[w2]], P[Sqrt[w1], Sqrt[w2]]}, {{w1, w2}}]

enter image description here

Or in the old coordinates, where I assume that r1>0 and r2>0:

Simplify[D[{Q[Sqrt[w1], Sqrt[w2]], 
    P[Sqrt[w1], Sqrt[w2]]}, {{w1, w2}}] /. {w1 -> r1^2, 
   w2 -> r2^2}, {r1 > 0, r2 > 0}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks for the insight! So you suggestion is to just use your last piece of code to compute the full Jacobian J? $\endgroup$
    – ofir_13
    Commented May 26, 2021 at 9:12
  • $\begingroup$ I think I might have solved the problem using the "substitution method" you showed above, very nice. I basically replaced the last line of the code I posted in the question with ``` FullSimplify[J[RTm, Rif] /. {RTm -> ((230 1.73)^4 + 4 (230 1.73)^2 0.2 100 3.14 Tm)/(4 (0.2)^2), Rif -> ((230 1.73) mif 100 3.14 )/((0.2)^2 + (100 3.14 0.0072)^2)}] ``` $\endgroup$
    – ofir_13
    Commented May 26, 2021 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.