Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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>