0
$\begingroup$

I have a very long code that take an extremely large amount of time. I have been using the While loop, but with each iteration, it gets slower and slower, so I know there has to be a better way of doing this than using While. I won't post the code since it's too long, but I will post the basic idea of what I need it to do.

currentvalue = 50; 
a = 5; 
b = 10; 
table1 = {10, 10, 5, 10, 5}; 
table2 = {10, 10, 5, 10, 5}; 
n = 1; While[n <= 15, rnd = RandomChoice[{1, 2, 3, 4, 5}]; If[table1[[rnd]] == a, table1[[rnd]] = b, table1[[rnd]] = a]; If[Total[table1] < currentvalue, currentvalue = Total[table1]; table2 = table1, table1 = table2]; Print[{table1, currentvalue}]; n++]

So, in this example, I am trying to minimize the sum of the list elements, so I am flipping one element to the other value, and checking if the sum gets smaller, if so, I update table1, and so on. So the idea is, I need to do many iterations, where I use the first conditional to initially update the list table1, and then I use the second conditional to see if I want to accept that update or not. Is there a way I can use Table with something like this so I can have my code run much faster than it is?

Edit: I will include a shortened version of the actual code in case it helps.

h = 1/10^2;
p = 1600;
stiffness = {2*10^6, 2*10^6, 2*10^6, 2*10^6};
boundary1 = .25;
boundary2 = boundary1*2;
boundary3 = boundary1*3;
boundary4 = boundary1*4;
a = 8*10^4;
b = 2*10^6;

function := 
 Module[{nsolution1, radialstretch1, thetastretch1, angle1, 
   radialboundary1, thetastretchderivative1, Rderivative1, f1, 
   nsolution2, radialstretch2, thetastretch2, angle2, radialboundary2,
    thetastretchderivative2, Rderivative2, f2, nsolution3, 
   radialstretch3, thetastretch3, angle3, radialboundary3, 
   thetastretchderivative3, Rderivative3, f3, nsolution4, 
   radialstretch4, thetastretch4, angle4, radialboundary4, 
   thetastretchderivative4, Rderivative4, f4, slope1, slope2, slope3, 
   slope4, , values, slopeerror, w1}, 
  nsolution1 = 
   ParametricNDSolve[{lt[r] + r Derivative[1][lt][r] == 
      lr[r] Cos[\[Beta][r]], 
     p r lt[r]* lr[r] Cos[\[Beta][r]] - 
       stiffness[[1]]* 
        h Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[1]]* 
        h r Derivative[1][\[Beta]][
         r] Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[1]]* 
        h r Sin[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) == 
      0, -p r lr[r] lt[r] Sin[\[Beta][r]] - 
       stiffness[[1]]* 
        h Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) + 
       stiffness[[1]]* 
        h r Derivative[1][\[Beta]][
         r] Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[1]]* 
        h r Cos[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) + 
       stiffness[[1]]* h (lt[r] - 1/(lr[r]^2 lt[r]^3)) == 0, 
     lt[0.001] == w1, lr[0.001] == w1, \[Beta][0.001] == 0.001}, {lr, 
     lt, \[Beta]}, {r, 0.001, boundary1}, {w1}];
  radialstretch1[v1_?NumericQ, v2_?NumericQ] := 
   lr[v1][v2] /. nsolution1; 
  thetastretch1[v1_?NumericQ, v2_?NumericQ] = 
   lt[v1][v2] /. nsolution1; 
  angle1[v1_?NumericQ, v2_?NumericQ] = \[Beta][v1][v2] /. nsolution1;
  radialboundary1[v1_?NumericQ, 
    v2_?NumericQ] := (1/(\[Sqrt]2)) (\[Sqrt]((stiffness[[
            1]]*(radialstretch1[v1, v2]^2 - 
             1/(radialstretch1[v1, v2]^2 thetastretch1[v1, v2]^2)))/
         stiffness[[
          2]] + (\[Sqrt]((stiffness[[1]] - 
               radialstretch1[v1, v2]^4 thetastretch1[v1, 
                 v2]^2 stiffness[[1]])^2 + 
             4 radialstretch1[v1, v2]^4 thetastretch1[v1, 
               v2]^2 (stiffness[[2]])^2))/(radialstretch1[v1, 
            v2]^2 thetastretch1[v1, v2]^2 stiffness[[2]])));
  thetastretchderivative1[v1_?NumericQ, v2_?NumericQ] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((\(lt[v1]\)[t] /.  
       nsolution1)\)\) /. t -> v2;
  Rderivative1[v1_?NumericQ, v2_?NumericQ] := 
   thetastretch1[v1, v2] + v2* thetastretchderivative1[v1, v2];
  f1[v1_?NumericQ, 
    v2_?NumericQ] := -\[Sqrt](radialstretch1[v1, v2]^2 - 
       Rderivative1[v1, v2]^2);
  slope1[v1_?NumericQ, v2_?NumericQ] := 
   1/Rderivative1[v1, v2]/f1[v1, v2];
  
  nsolution2 = 
   ParametricNDSolve[{lt[r] + r Derivative[1][lt][r] == 
      lr[r] Cos[\[Beta][r]], 
     p r lt[r] lr[r] Cos[\[Beta][r]] - 
       stiffness[[2]]* 
        h Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[2]]* 
        h r Derivative[1][\[Beta]][
         r] Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[2]]* 
        h r Sin[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) == 
      0, -p r lr[r] lt[r] Sin[\[Beta][r]] - 
       stiffness[[2]]* 
        h Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) + 
       stiffness[[2]]* 
        h r Derivative[1][\[Beta]][
         r] Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[2]]* 
        h r Cos[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) + 
       stiffness[[2]]* h (lt[r] - 1/(lr[r]^2 lt[r]^3)) == 0, 
     lr[boundary1] == radialboundary1[w1, boundary1], 
     lt[boundary1] == 
      thetastretch1[w1, boundary1], \[Beta][boundary1] == 
      angle1[w1, boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, 
     boundary2}, {w1}];
  radialstretch2[v1_?NumericQ, v2_?NumericQ] := 
   lr[v1][v2] /. nsolution2; 
  thetastretch2[v1_?NumericQ, v2_?NumericQ] = 
   lt[v1][v2] /. nsolution2; 
  angle2[v1_?NumericQ, v2_?NumericQ] = \[Beta][v1][v2] /. nsolution2;
  radialboundary2[v1_?NumericQ, 
    v2_?NumericQ] := (1/(\[Sqrt]2)) (\[Sqrt]((stiffness[[
            2]]*(radialstretch2[v1, v2]^2 - 
             1/(radialstretch2[v1, v2]^2 thetastretch2[v1, v2]^2)))/
         stiffness[[
          3]] + (\[Sqrt]((stiffness[[2]] - 
               radialstretch2[v1, v2]^4 thetastretch2[v1, 
                 v2]^2 stiffness[[2]])^2 + 
             4 radialstretch2[v1, v2]^4 thetastretch2[v1, 
               v2]^2 (stiffness[[3]])^2))/(radialstretch2[v1, 
            v2]^2 thetastretch2[v1, v2]^2 stiffness[[3]])));
  thetastretchderivative2[v1_?NumericQ, v2_?NumericQ] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((\(lt[v1]\)[t] /.  
       nsolution2)\)\) /. t -> v2;
  Rderivative2[v1_?NumericQ, v2_?NumericQ] := 
   thetastretch2[v1, v2] + v2* thetastretchderivative2[v1, v2];
  f2[v1_?NumericQ, 
    v2_?NumericQ] := -\[Sqrt](radialstretch2[v1, v2]^2 - 
       Rderivative2[v1, v2]^2);
  slope2[v1_?NumericQ, v2_?NumericQ] := 
   1/Rderivative2[v1, v2]/f2[v1, v2];
  nsolution3 = 
   ParametricNDSolve[{lt[r] + r Derivative[1][lt][r] == 
      lr[r] Cos[\[Beta][r]], 
     p r lt[r] lr[r] Cos[\[Beta][r]] - 
       stiffness[[3]]* 
        h Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[3]]* 
        h r Derivative[1][\[Beta]][
         r] Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[3]]* 
        h r Sin[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) == 
      0, -p r lr[r] lt[r] Sin[\[Beta][r]] - 
       stiffness[[3]]* 
        h Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) + 
       stiffness[[3]]* 
        h r Derivative[1][\[Beta]][
         r] Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[3]]* 
        h r Cos[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) + 
       stiffness[[3]]* h (lt[r] - 1/(lr[r]^2 lt[r]^3)) == 0, 
     lr[boundary2] == radialboundary2[w1, boundary2], 
     lt[boundary2] == 
      thetastretch2[w1, boundary2], \[Beta][boundary2] == 
      angle2[w1, boundary2]}, {lr, lt, \[Beta]}, {r, boundary2, 
     boundary3}, {w1}];
  radialstretch3[v1_?NumericQ, v2_?NumericQ] := 
   lr[v1][v2] /. nsolution3; 
  thetastretch3[v1_?NumericQ, v2_?NumericQ] = 
   lt[v1][v2] /. nsolution3; 
  angle3[v1_?NumericQ, v2_?NumericQ] = \[Beta][v1][v2] /. nsolution3;
  radialboundary3[v1_?NumericQ, 
    v2_?NumericQ] := (1/(\[Sqrt]2)) (\[Sqrt]((stiffness[[
            3]]*(radialstretch3[v1, v2]^2 - 
             1/(radialstretch3[v1, v2]^2 thetastretch3[v1, v2]^2)))/
         stiffness[[
          4]] + (\[Sqrt]((stiffness[[3]] - 
               radialstretch3[v1, v2]^4 thetastretch3[v1, 
                 v2]^2 stiffness[[3]])^2 + 
             4 radialstretch3[v1, v2]^4 thetastretch3[v1, 
               v2]^2 (stiffness[[4]])^2))/(radialstretch3[v1, 
            v2]^2 thetastretch3[v1, v2]^2 stiffness[[4]])));
  thetastretchderivative3[v1_?NumericQ, v2_?NumericQ] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((\(lt[v1]\)[t] /.  
       nsolution3)\)\) /. t -> v2;
  Rderivative3[v1_?NumericQ, v2_?NumericQ] := 
   thetastretch3[v1, v2] + v2* thetastretchderivative3[v1, v2];
  f3[v1_?NumericQ, 
    v2_?NumericQ] := -\[Sqrt](radialstretch3[v1, v2]^2 - 
       Rderivative3[v1, v2]^2);
  slope3[v1_?NumericQ, v2_?NumericQ] := 
   1/Rderivative3[v1, v2]/f3[v1, v2];
  nsolution4 = 
   ParametricNDSolve[{lt[r] + r Derivative[1][lt][r] == 
      lr[r] Cos[\[Beta][r]], 
     p r lt[r] lr[r] Cos[\[Beta][r]] - 
       stiffness[[4]]* 
        h Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[4]]* 
        h r Derivative[1][\[Beta]][
         r] Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[4]]* 
        h r Sin[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) == 
      0, -p r lr[r] lt[r] Sin[\[Beta][r]] - 
       stiffness[[4]]* 
        h Cos[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) + 
       stiffness[[4]]* 
        h r Derivative[1][\[Beta]][
         r] Sin[\[Beta][r]] (lr[r] - 1/(lr[r]^3 lt[r]^2)) - 
       stiffness[[4]]* 
        h r Cos[\[Beta][
          r]] (Derivative[1][lr][
           r] + (3 Derivative[1][lr][r])/(lr[r]^4 lt[
              r]^2) + (2 Derivative[1][lt][r])/(lr[r]^3 lt[r]^3)) + 
       stiffness[[4]]* h (lt[r] - 1/(lr[r]^2 lt[r]^3)) == 0, 
     lr[boundary3] == radialboundary3[w1, boundary3], 
     lt[boundary3] == 
      thetastretch3[w1, boundary3], \[Beta][boundary3] == 
      angle3[w1, boundary3]}, {lr, lt, \[Beta]}, {r, boundary3, 
     boundary4}, {w1}];
  radialstretch4[v1_?NumericQ, v2_?NumericQ] := 
   lr[v1][v2] /. nsolution4; 
  thetastretch4[v1_?NumericQ, v2_?NumericQ] = 
   lt[v1][v2] /. nsolution4; 
  angle4[v1_?NumericQ, v2_?NumericQ] = \[Beta][v1][v2] /. nsolution4;
  thetastretchderivative4[v1_?NumericQ, v2_?NumericQ] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((\(lt[v1]\)[t] /.  
       nsolution4)\)\) /. t -> v2;
  Rderivative4[v1_?NumericQ, v2_?NumericQ] := 
   thetastretch4[v1, v2] + v2* thetastretchderivative4[v1, v2];
  f4[v1_?NumericQ, 
    v2_?NumericQ] := -\[Sqrt](radialstretch4[v1, v2]^2 - 
       Rderivative4[v1, v2]^2);
  slope4[v1_?NumericQ, v2_?NumericQ] := 
   1/Rderivative4[v1, v2]/f4[v1, v2];
  values = Table[x, {x, 0.05, 1, 0.05}];
  slopeerror[
    v1_?NumericQ] := (1/
      Length[values])*(Sum[(slope1[v1, values[[i]]] + 
         2*values[[i]]*Exp[-values[[i]]^2])^2, {i, 1, 5}] + 
      Sum[(slope2[v1, values[[i]]] + 
         2*values[[i]]*Exp[-values[[i]]^2])^2, {i, 5, 10}] + 
      Sum[(slope3[v1, values[[i]]] + 
         2*values[[i]]*Exp[-values[[i]]^2])^2, {i, 10, 15}] + 
      Sum[(slope4[v1, values[[i]]] + 
         2*values[[i]]*Exp[-values[[i]]^2])^2, {i, 15, 20}]);
  
  
  {thetastretch4, slopeerror}]

current = 200;
stiffness0 = {8*10^4, 2*10^6, 8*10^4, 2*10^6, 8*10^4, 2*10^6, 8*10^4, 
   2*10^6, 8*10^4, 2*10^6, 8*10^4};

Quiet[Monitor[n = 1; 
  While[n < 25, 
   rnd = RandomChoice[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}]; 
   If[stiffness[[rnd]] == a, stiffness[[rnd]] = b, 
    stiffness[[rnd]] = a];
   listofvalues = 
    Quiet[Select[
      Table[{x, 
        function[[1]][x, 1]}, {x, {1.001, 1.002, 1.003, 1.004, 1.005, 
         1.0055, 1.006, 1.007, 1.008, 1.009, 1.01, 1.02, 1.03, 1.04, 
         1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 
         1.15, 1.16, 1.17, 1.18, 1.2, 1.22, 1.24, 1.26, 1.28, 1.3, 
         1.32, 1.34, 1.36, 1.38, 1.4}}], Abs[#[[2]]] < 5 &]];
   root = 
    w1 /. FindRoot[
      function[[1]][w1, 1] == 1, {w1, listofvalues[[1, 1]], 
       listofvalues[[3, 1]]}]; 
   If[function[[2]][root] < current, current = function[[2]][root]; 
    stiffness0 = stiffness; 
    Print[function[[2]][root], stiffness, root], 
    stiffness = stiffness0;]; n++]; n, n]]

So, I have these differential equations that I defined inside of Module in what I called function, and I need them solved for different conditions many times over. That is why I am using that while loop, where I change the list called stiffness, and then redo the equations to see if I am getting closer to a desired result. Each time the stiffness list is changed, "function" gives a different result for the differential equations, so it has to be ran again inside the loop.

Improve this question
$\endgroup$
7
  • $\begingroup$ I think you could get better answers if you would describe the problem you are attempting to solve because there is probably a significantly better abstraction available in Mathematica for it than While loops or brute-force Tables. I believe current code in the question doesn't sufficiently expose the nature of the problem to write informative answers. $\endgroup$
    – kirma
    Commented Mar 4, 2023 at 7:52
  • $\begingroup$ I will edit the question to include the code as well. $\endgroup$
    – juv95
    Commented Mar 4, 2023 at 7:58
  • $\begingroup$ Even better if you can describe the problem you're trying to solve; for instance "find a subset of size 5 with a minimal sum." If the problem itself is well-formed there tend to be good domain-specific tools available in Mathematica... $\endgroup$
    – kirma
    Commented Mar 4, 2023 at 8:04
  • $\begingroup$ I've updated the question with the edit and a short explanation of the problem itself. But, it is a bit hard to look at.. $\endgroup$
    – juv95
    Commented Mar 4, 2023 at 8:07
  • 1
    $\begingroup$ @DanielHuber I believe Total is a bit unfortunate choice for an arbitrary function in this case... $\endgroup$
    – kirma
    Commented Mar 4, 2023 at 17:57

1 Answer 1

Reset to default
1
$\begingroup$

Here's something resembling your problem in a functional programming form: a complicated function f which takes a vector as input, ab listing two possible values of each vector item, and a Nested search which flips individual vector items randomly, attempting to avoid repeated calculations of f. I can't say anything about efficiency, but I assume most of the time in the real-world case would be spent inside f anyway...

Module[
 {ab, f},
 (* Alternative values on each vector item. *)
 ab = {1, Sqrt[2]};
 (* A complicated function we want to minimize. *)
 f[v_] := Sum[Sin[Prime[i] v[[i]]], {i, Length[v]}] // N;
 Nest[
   (* Search iteration step. *)
   Apply[Function[{v, best, seen},
     With[
      {vcand =
        (* Find a candidate vector by flipping a random vector entry.
           If vector is already tested, try again until a new one
           is found, or Length[v] tries have been made. *)
        NestWhile[
         ReplaceAt[v, {1 -> 2, 2 -> 1}, RandomInteger[Length[v]]] &,
         {},
         # == {} || MemberQ[seen, #] &,
         1, Length[v]]},
      With[{val = f[ab[[vcand]]]},
       If[val <= best,
        (* If value was improved pass new candidate to
           the next iteration. *)
        {vcand, Echo[val], Append[seen, vcand]},
        {v, best, Append[seen, vcand]}]]]]],
   (* Start with random vector, infinite best value and
      no vectors seen.*)
   {RandomInteger[{1, 2}, {10}], Infinity, {}},
   (* 100 iterations *)
   100] //
  (* Assign ab values to the found result. *)
  ab[[First[#]]] &]

(* 
-2.38285

-4.1252

-5.27623

-5.76921

-5.76921

-5.76921

-5.76921

-5.76921

-5.76921

-5.76921

-5.76921

-5.76921
*)

(* {Sqrt[2], Sqrt[2], 1, Sqrt[2], 1, Sqrt[2], 1, 1, 1, 1} *)

The minimum in this case converged to the real minimum:

Sum[Sin[Prime[i] #[[i]]], {i, Length[#]}] & /@ 
   Tuples[{1, Sqrt[2]}, 10] // Min // N

(* -5.76921 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ I haven’t had the chance to try it in my problem, but it seems like this is exactly what I was asking for. Thank you! $\endgroup$
    – juv95
    Commented Mar 6, 2023 at 17:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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