3

I use the Runge Kutta method in Mathematica to solve a matrix differential equation, that depends on a number of variables. At the end of it, I extract some information that depends on my choice of variables at the start.

I want to be able to numerically optimise this such that the information I extract at the end is done with the best possible choice of variables. A fairly basic version of what I want to accomplish is written below

    opt[var_] := (
    cc1 =({{var/2, 1},{0, 1/(2 var)}});
    dd1 = ( {{0, t/var},{-var t, t}});

    rk={};
    it = 0; 
    steps = 10^4; 
    dt = 10/steps; 

    sf =cc1;

    lya[mat_, time_] := (dd1 /. t -> time).mat + mat.(dd1 /. t -> time)\[Transpose];

      Do[si = sf // N;
         k1 = lya[si, it] // N;
         k2 = lya[si + dt/2 k1, it + dt/2] // N;
         k3 = lya[si + dt/2 k2, it + dt/2] // N;
         k4 = lya[si + dt k3, it] // N;

       sf = si + \dt/6 (k1 + 2 k2 + 2 k3 + k4) // N;
      AppendTo[rk, {it, sf}];, {iterations, steps}];

     ee = var/2 sf[[1, 1]] + sf[[2, 2]]/(2 var)
    )

Let's say ee is what I want to optimise, and var is my variable. I assumed I could just call

    NMinimize[{opt[var],0<var<=5},var]

but the RungeKutta won't solve. I've looked around the Mathematica documentation for something around this, and came up with naught.

Any help is appreciated. More information available if required.

For context the actual code I'll be running has 6x6 matrices, with 10 variables to optimise over! these may be able to be reduced.

Thanks

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Brendan Reid
  • 58
  • 7
  • 6
    Is a self implementation of RK4 necessary? Why not use e.g. ParametricNDSolveValue to solve the equation? – xzczd Feb 18 '17 at 02:33
  • 2
    The code contains \dt/6, which is a syntax error Also, it is initialized as 0 but never changes, even though it is used throughout the remainder of the code and saved in rk. This, too, probably is an error. Please correct. – bbgodfrey Feb 18 '17 at 05:09
  • dt in my original code is actually delta, thus the erroneous \ when I was editing the post. & youre correct, it is initialised as 0 and never changes in this code. This particular block was written solely for this post, as an example of what I need it to do. There should be an extra it = it+ dt inside the Do. – Brendan Reid Feb 18 '17 at 11:04
  • You can click the edit button to correct the mistake in your post. – xzczd Feb 22 '17 at 02:15

1 Answers1

1

You need to define opt with NumericQ,

opt[var_?NumericQ] := ...

to keep the function from being evaluated symbolically, which does not finish for a long, long time. Your opt function works for me on numeric inputs, but your routine takes a long time:

opt[2.] // AbsoluteTiming
(* {1.52733, 1.06291} *)

That's just for one evaluation. I didn't time this, but let it run over dinner:

FindMinimum[{opt[var], 0 < var <= 5}, {var, 2.}]
(* {0.501664, {var -> 1.00167}} *)

NMinimize[] would work, too, but will take much longer, most likely, given that it tends to evaluate the objective function more times.

Community
  • 1
Goofy
  • 2,767
  • 10
  • 12