Re: Problem evaluating the Jacobian
- To: mathgroup at smc.vnet.net
- Subject: [mg123886] Re: Problem evaluating the Jacobian
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Wed, 28 Dec 2011 05:15:59 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201112262230.RAA20917@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
> f[x_]:={x[[1]]...........};
> df[x_]:={{x[[1]]...},{....},.....};
>
> The variables have to be written as part of a list.
Those are not variables. They're the parts of a list.
If you want a three-dimensional list of variables, you might try
Clear[x]
Array[x, {3, 2, 4}]
{{{x[1, 1, 1], x[1, 1, 2], x[1, 1, 3], x[1, 1, 4]}, {x[1, 2, 1],
x[1, 2, 2], x[1, 2, 3], x[1, 2, 4]}}, {{x[2, 1, 1], x[2, 1, 2],
x[2, 1, 3], x[2, 1, 4]}, {x[2, 2, 1], x[2, 2, 2], x[2, 2, 3],
x[2, 2, 4]}}, {{x[3, 1, 1], x[3, 1, 2], x[3, 1, 3],
x[3, 1, 4]}, {x[3, 2, 1], x[3, 2, 2], x[3, 2, 3], x[3, 2, 4]}}}
Bobby
On Mon, 26 Dec 2011 16:30:25 -0600, Joao <joaopereira9 at gmail.com> wrote:
> Hi there,
>
> I need some help with the following matter:
>
> I am in the middle of building an implicit runge-kutta method (stage 3),
> and at sometime in the process, a newton algorithm will have to be used.
> I have built a customized newton algorithm that takes the function and
> jacobian as arguments, but they have to be in the format:
>
> f[x_]:={x[[1]]...........};
> df[x_]:={{x[[1]]...},{....},.....};
>
> The variables have to be written as part of a list.
>
> This format restriction and the need to have the implicit runge - kutta
> method flexible enough to be used in any dimension, has led me to
> construct a few intermediate steps in order to have the function and the
> jacobian to be used in the newton algorithm (inside the runge-kutta
> algorithm) working properly. Note that if the original function has n
> dimensions, the function to be used in the newton algorithm will have
> 3*n dimensions. The unfinished runge-kutta algorithm is
>
> irk[frk_, t0_, tf_, x0_, irkn_] :=
> Module[{t0i, x0i, irkh, a, b, c, ni, t1, t2, t3, uarg, u1, u2, u3,
> u123, fi, var, yeq, dyeq, v},
> Clear[t0i, x0i, yeq, v, var, dyeq];
> t0i = N[t0];
> x0i = N[x0];
> irkh = (tf - t0)/(irkn -
> 1);(*Método implicito de 3 estágios com ordem 5, Radau IIA*)
> a = N[{{(88. - 7*6^0.5)/360, (296. - 169*6^0.5)/
> 1800, (-2. + 3*6^0.5)/225}, {(296. + 169*6^0.5)/
> 1800, (88. + 7*6^0.5)/360, (-2. - 3*6^0.5)/
> 225}, {(16. - 6^0.5)/36, (16. + 6^0.5)/36, 1/9}}];
> b = N[{(16. - 6.^0.5)/36, (16. - 6.^0.5)/36, 1/9}];
> c = N[{(4. - 6.^0.5)/10, (4. + 6.^0.5)/10, 1}];
> ni = Length[Apply[frk, {t0, x0}]];
> fi[v_] = Array[v[[#2 + (#1 - 1)*ni]] &, {3, ni}];
> var[v_] =
> Join[Transpose[{fi[v][[1]]}], Transpose[{fi[v][[2]]}],
> Transpose[{fi[v][[3]]}]][[All, 1]];
> t1 = t0i + c[[1]]*irkh;
> t2 = t0i + c[[2]]*irkh;
> t3 = t0i + c[[3]]*irkh;
> uarg[v_] =
> Join[Transpose[{(a.fi[v])[[1]]}], Transpose[{(a.fi[v])[[2]]}],
> Transpose[{(a.fi[v])[[3]]}]][[All, 1]];
> u1[v_] = x0i + irkh*(uarg[v][[1 ;; ni]]);
> u2[v_] = x0i + irkh*(uarg[v][[ni + 1 ;; 2 ni]]);
> u3[v_] = x0i + irkh*(uarg[v][[2 ni + 1 ;; 3 ni]]);
> u123[v_] =
> Join[Apply[frk, {t1, u1[v]}], Apply[frk, {t2, u2[v]}],
> Apply[frk, {t3, u3[v]}]]; Print["u123 ", u123[v]];
> yeq[v_] := {var[v] - u123[v]}[[1]];
> Print["yeq[v_] ", MatrixForm[yeq[v]]];
> dyeq[v_] := Outer[D, yeq[v], var[v]];
> Print["dyeq[v] ", MatrixForm[dyeq[v]]];
> Print["function ", Apply[yeq, {{1., 2., 0.5, 5., 2.5, 0.6}}]];
> Print["Jacobian ", Apply[dyeq, {{1., 2., 0.5, 5., 2.5, 0.6}}]];
> ];
>
> Now, it seems it is doing everything I want to up to the evaluation of
> the jacobian (arbitrary values). The evaluation of the Jacobian has to
> be done with the function Apply as this is the one used in the newton
> algorithm, and I do not want to change that algorithm as it is already
> working correctly.
>
> As a very simple example, I am using the function describing the right
> hand side of the differential equations
>
> f[t_, y_] := {y[[1]]^2 - y[[2]], y[[1]] + 0.01 y[[2]]};
>
> Any ideas why this is not working?
>
>
> Thanks in advance
>
> Joao Pereira
>
--
DrMajorBob at yahoo.com
- References:
- Problem evaluating the Jacobian
- From: Joao <joaopereira9@gmail.com>
- Problem evaluating the Jacobian