2
$\begingroup$

At the moment I have implemented the code for a Taylor 2nd order series for the function in three variables:

$x_3^3+\frac{x_1-x_2}{x_1+x_2}$

The code builds on following expression:

enter image description here

ClearAll["Global`*"]
Remove["Global`*"]

thetacalc[xtem_, utem_, nn_] := {

n = nn;
xx = Table[Subscript[x, i], {i, 1, n}];
uu = Table[Subscript[u, i], {i, 1, n}];

f1[x_] := ((x[[1]] - x[[2]])/(x[[1]] + x[[2]])) + (x[[3]])^3;
hh = Simplify[D[f1[xx + θ*uu], {xx, 2}]] // MatrixForm;

f2[f_, x_, u_] := 
f + Sum[uu[[j]]*D[f, {xx[[j]]}], {j, 1, n}] + 
1/2*Sum[Sum[uu[[i]]*uu[[j]]*hh[[1, i, j]], {j, 1, n}], {i, 1, n}];

sol2 = FullSimplify[f2[f1[xx], xx, uu]];
uu = utem;
xx = xtem;
subst = 
Flatten[{Table[Subscript[x, i] -> xx[[i]], {i, 1, n}], 
Table[Subscript[u, i] -> uu[[i]], {i, 1, n}]}];
subst2 = Flatten[{Table[Subscript[x, i] -> xx[[i]], {i, 1, n}]}];
diff = f1[xx + uu] - sol2 /. subst;
{sei = NSolve[diff == 0, θ, Reals], diff, FullSimplify[hh] /. subst2, sein = FullSimplify[hh] /. Flatten[{subst, sei}], Eigenvalues[sein[[1]]]}

}

thetacalc[{10, 5, 1}, {1, 0, 0}, 3] 

The output (theta, difference-function, H-Matrix [unevaluated], [evaluated], Eigenvalues) for x and u as given in thetacalc:

$\tiny \begin{array}{ccccc} \{\{\theta \to 0.326189\}\} & \frac{10}{(\theta +15)^3}-\frac{1}{360} & \left( \begin{array}{ccc} -\frac{4 \left(\theta u_2+5\right)}{\left(\theta \left(u_1+u_2\right)+15\right){}^3} & \frac{2 \left(\theta u_1-\theta u_2+5\right)}{\left(\theta \left(u_1+u_2\right)+15\right){}^3} & 0 \\ \frac{2 \left(\theta u_1-\theta u_2+5\right)}{\left(\theta \left(u_1+u_2\right)+15\right){}^3} & \frac{4 \left(\theta u_1+10\right)}{\left(\theta \left(u_1+u_2\right)+15\right){}^3} & 0 \\ 0 & 0 & 6 \left(\theta u_3+1\right) \\ \end{array} \right) & \left( \begin{array}{ccc} -0.00555556 & 0.00295899 & 0 \\ 0.00295899 & 0.0114735 & 0 \\ 0 & 0 & 6. \\ \end{array} \right) & \{6.,0.011973,-0.00605506\} \\ \end{array}$

The remainder will be adjusted with a theta between (0,1) from the mean value theorem to equal the original function (at point x with direction vector). A theta is given if the direction vector u is unequal to null-vector.

I am wondering if one could do this simpler as to avoid manual term calculation like the Hesse-matrix (in my code: hh) or even multiple sums.

Thank you for your time.

$\endgroup$
3
  • $\begingroup$ A side note, try avoiding the use of Subscript for variables, when you do Subscript[x, i]=1 you are actually saying Set a downvalue to Subscript not to x. $\endgroup$
    – rhermans
    Jan 12, 2016 at 17:06
  • $\begingroup$ Is is then not at all recommended to use indexed variables like above? For visibility I would like to still use them for not having to use every letter on my keyboard. Could you give me an alternative? Thanks $\endgroup$ Jan 12, 2016 at 19:13
  • $\begingroup$ Have a read in this question. $\endgroup$
    – rhermans
    Jan 12, 2016 at 19:31

1 Answer 1

3
$\begingroup$

Probably:

p = {x, y, z};
ff[x_, y_, z_] := (x - y)/(x + y) + z^3
hhh = D[ff[x, y, z], {p, 2}] // Simplify
h1 = hhh /. MapThread[Rule, {p, p + t { u1, u2, u3}}] /. 
               Thread[p :> {10, 5, 1}]

Mathematica graphics

You can get the same result you got with your code with the following (please note that most of the code is just formatting,not sure why you may want that on a function):

calc[ff_, p0_, u0_] := Module[{p = {x, y, z}, u = {u1, u2, u3}, c1, c2, t, diff, tval, mt},
  mt = Flatten@MapThread[Rule, {{p, u}, {p0, u0}}, 2];
  c1 = D[ff@p, {p, 1}] // Simplify;
  c2 = D[ff[p + t u], {p, 2}] // Simplify;
  diff = (ff[p + u] - ff[p] - u.c1 - 1/2 u.c2.u) /. mt;
  tval = First@NSolve[diff == 0, Reals];

  (*Output Formatting follows*)
  {tval /. t -> \[FormalT],
   diff /. t -> \[FormalT],
   c2 /. Thread[p :> p0] /. t -> \[FormalT] // MatrixForm,
   c2 /. mt /. tval // MatrixForm,
   Eigenvalues[c2] /. mt /. tval}]

Usage

f[{x_, y_, z_}] := (x - y)/(x + y) + z^3
calc[f, {10, 5, 1}, {1, 0, 0}]

Mathematica graphics

$\endgroup$
2
  • $\begingroup$ You may want to use formal symbols for {x, y, z, u1,u2, u3} to avoid variable scoping issues. Or use Unique[ ] if you don't need to display the unevaluated forms nicely $\endgroup$ Jan 13, 2016 at 6:59
  • 1
    $\begingroup$ This is great work for getting a more concise code/module. A higher order generalization would be easier to implement in here as well. Building on this module on could also extend it for optimization. $\endgroup$ Jan 14, 2016 at 15:09

Your Answer

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

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