Services & ResourcesWolfram ForumsMathGroup Archive

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

Problem evaluating the Jacobian

• To: mathgroup at smc.vnet.net
• Subject: [mg123879] Problem evaluating the Jacobian
• From: Joao <joaopereira9 at gmail.com>
• Date: Mon, 26 Dec 2011 17:30:25 -0500 (EST)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com

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

• Follow-Ups:
  • Re: Problem evaluating the Jacobian
    • From: DrMajorBob <btreat1@austin.rr.com>