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$
  • $\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 '16 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$ – chrisoutwright Jan 12 '16 at 19:13
  • $\begingroup$ Have a read in this question. $\endgroup$ – rhermans Jan 12 '16 at 19:31
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

| improve this answer | |
$\endgroup$
  • $\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$ – Dr. belisarius Jan 13 '16 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$ – chrisoutwright Jan 14 '16 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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