# Generalize the code to more variables

I have this code it runs and gives me the solution. How can I make it more compact and If I want it to extend(generalize) it to more variables how can I do it. I had asked a similar question here (Extend the code for more variables (at least 100)) but I am unable to apply that to these equations. The first problem is to create equations for any arbitrary variable. I want to extend this to around 50 variables I can type the equations but I am assuming that there will be a better way to generalize. Any help is highly appreciated.

It has $$x0=0$$ and $$x10=1$$ and 9 unknowns.

    f1[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x1) -
100/201 ((x2) - 1/100 Cos[(5 (x2))^2])
f2[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x2) -
100/201 ((x1) + (x3) - 1/100 Cos[(5 (x3) - 5 (x1))^2])
f3[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x3) -
100/201 ((x2) + (x4) - 1/100 Cos[(5 (x4) - 5 (x2))^2])
f4[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x4) -
100/201 ((x3) + (x5) - 1/100 Cos[(5 (x5) - 5 (x3))^2])
f5[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x5) -
100/201 ((x4) + (x6) - 1/100 Cos[(5 (x6) - 5 (x4))^2])
f6[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x6) -
100/201 ((x5) + (x7) - 1/100 Cos[(5 (x7) - 5 (x5))^2])
f7[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x7) -
100/201 ((x6) + (x8) - 1/100 Cos[(5 (x8) - 5 (x6))^2])
f8[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x8) -
100/201 ((x7) + (x9) - 1/100 Cos[(5 (x9) - 5 (x7))^2])
f9[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] := (x9) -
100/201 ((x8) + 1 - 1/100 Cos[(5 - 5 (x8))^2])

f[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_,
x9_] := {f1[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f2[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f3[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f4[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f5[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f6[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f7[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f8[x1, x2, x3, x4, x5, x6, x7, x8, x9],
f9[x1, x2, x3, x4, x5, x6, x7, x8, x9]}
g1[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x2) - 1/100 Cos[(5 (x2))^2])
g2[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x1) + (x3) - 1/100 Cos[(5 (x3) - 5 (x1))^2])
g3[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x2) + (x4) - 1/100 Cos[(5 (x4) - 5 (x2))^2])
g4[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x3) + (x5) - 1/100 Cos[(5 (x5) - 5 (x3))^2])
g5[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x4) + (x6) - 1/100 Cos[(5 (x6) - 5 (x4))^2])
g6[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x5) + (x7) - 1/100 Cos[(5 (x7) - 5 (x5))^2])
g7[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x6) + (x8) - 1/100 Cos[(5 (x8) - 5 (x6))^2])
g8[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x7) + (x9) - 1/100 Cos[(5 (x9) - 5 (x7))^2])
g9[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_] :=
100/201 ((x8) + 1 - 1/100 Cos[(5 - 5 (x8))^2])

g[x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_,
x9_] := {g1[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g2[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g3[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g4[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g5[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g6[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g7[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g8[x1, x2, x3, x4, x5, x6, x7, x8, x9],
g9[x1, x2, x3, x4, x5, x6, x7, x8, x9]}

jac[ {x1, x2, x3, x4, x5, x6, x7, x8, x9}] = {D[
g1[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g2[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g3[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g4[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g5[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g6[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g7[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g8[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}],
D[g9[x1, x2, x3, x4, x5, x6, x7, x8,
x9], {{x1, x2, x3, x4, x5, x6, x7, x8, x9}}]};

{t1, t2, t3, t4, t5, t6, t7, t8, t9} = {1/5, 1/5, 1/5, 1/5, 1/5, 1/5,
1/5, 1/5, 1/5};
eps = 10^-100;
max = 10;
per = 200;

a = 1;
b = SetAccuracy[Norm[f[t1, t2, t3, t4, t5, t6, t7, t8, t9], Infinity],
per];
k = 1;
While[k <= max && a + b > eps,
M = SetAccuracy[-
jac[ {x1, x2, x3, x4, x5, x6, x7, x8, x9}] /. {x1 -> t1,
x2 -> t2, x3 -> t3, x4 -> t4, x5 -> t5, x6 -> t6, x7 -> t7,
x8 -> t8, x9 -> t9}, per];
J = SetAccuracy[
Inverse[IdentityMatrix[9] -
jac[ {x1, x2, x3, x4, x5, x6, x7, x8, x9}] /. {x1 -> t1,
x2 -> t2, x3 -> t3, x4 -> t4, x5 -> t5, x6 -> t6, x7 -> t7,
x8 -> t8, x9 -> t9}], per];
{z1, z2, z3, z4, z5, z6, z7, z8, z9} =
SetAccuracy[
J . M . {t1, t2, t3, t4, t5, t6, t7, t8, t9} +
J . g[t1, t2, t3, t4, t5, t6, t7, t8, t9], per];
{y1, y2, y3, y4, y5, y6, y7, y8, y9} =
SetAccuracy[g[z1, z2, z3, z4, z5, z6, z7, z8, z9], per];
{h1, h2, h3, h4, h5, h6, h7, h8, h9} =
SetAccuracy[g[y1, y2, y3, y4, y5, y6, y7, y8, y9], per];
a = SetAccuracy[
Norm[({t1, t2, t3, t4, t5, t6, t7, t8, t9} - {h1, h2, h3, h4, h5,
h6, h7, h8, h9}), Infinity], per];
b = SetAccuracy[
Norm[f[h1, h2, h3, h4, h5, h6, h7, h8, h9], Infinity], per];
{t1, t2, t3, t4, t5, t6, t7, t8, t9} = {h1, h2, h3, h4, h5, h6, h7,
h8, h9};
Print["k=", NumberForm[k, 20]];
Print["{h1,h2,h3,h4,h5,h6,h7,h8,h9}=",
NumberForm[{h1, h2, h3, h4, h5, h6, h7, h8, h9}, 20]];
Print["a=", NumberForm[a, 20]];
Print["b=", NumberForm[b, 20]];
k = k + 1;];
Print["The number of iterations is:" NumberForm[(k - 1), 2]];
Print["a+b=" NumberForm[a + b, 20]];

• You could simplify the variable number of arguments by passing a list rather than individual arguments. i.e., for your first definition f1[parms_] := parms[[1]] - 100/201 (parms[[2]] - 1/100 Cos[(5*parms[[2]])^2]); f1[{1, 2, 3}]. You pass the list but the function only uses what it needs. Commented Oct 6, 2023 at 17:10
• It doesn't work Commented Oct 6, 2023 at 17:20

Redefine f and g to take a list of variables

Clear["Global*"]

f[vars_] := Module[{t1, t2}, Table[
If[n == 1, t1 = 0, t1 = vars[[n - 1]]];
If[n == Length[vars], t2 = 1, t2 = vars[[n + 1]]];
vars[[n]] - 100/201 (t1 + t2 - 1/100 Cos[(5 t2 - 5 t1)^2])
, {n, Length[vars]}]
]

g[vars_] := Module[{t1, t2}, Table[
If[n == 1, t1 = 0, t1 = vars[[n - 1]]];
If[n == Length[vars], t2 = 1, t2 = vars[[n + 1]]];
100/201 (t1 + t2 - 1/100 Cos[(5 t2 - 5 t1)^2])
, {n, Length[vars]}]
]

eps = 10^-100;
max = 10;
per = 200;
a = 1;

vars = Array[x, 9];
varst = Table[1/5, {Length[vars]}];
r = Table[vars[[n]] -> varst[[n]], {n, Length[vars]}];
jac2 = Table[D[gg, {vars}], {gg, g[vars]}];
k = 1;
b = SetAccuracy[Norm[f[varst], Infinity], per]

While[k <= max && a + b > eps,
M = SetAccuracy[-jac2 /. r, per];
J = SetAccuracy[Inverse[IdentityMatrix[9] - jac2 /. r], per];
varsz = SetAccuracy[J . M . varst + J . g[varst], per];
varsy = SetAccuracy[g[varsz], per];
varsh = SetAccuracy[g[varsy], per];
a = SetAccuracy[Norm[(varst - varsh), Infinity], per];
b = SetAccuracy[Norm[f[varsh], Infinity], per];
varst = varsh;
Print["k=", NumberForm[k, 20]];
Print["{h1,h2,h3,h4,h5,h6,h7,h8,h9}=", NumberForm[varsh, 20]];
Print["a=", NumberForm[a, 20]];
Print["b=", NumberForm[b, 20]];
k = k + 1;];
Print["The number of iterations is:" NumberForm[(k - 1), 2]];
Print["a+b=" NumberForm[a + b, 20]];


• Hi, Thanks for the beautiful code. However its not wrorking. The command jac2 is not working and I dont know that 'r' is when you do -jac2 /. r Commented Oct 7, 2023 at 7:24
• @Learner - Apologies, r is the list of replacement rules, which I forgot to paste. Answer updated to include r = Table[vars[[n]] -> varst[[n]], {n, Length[vars]}]; Commented Oct 7, 2023 at 17:45

Here's a start. Encode the general pattern:

newF[a_, b_, c_] := b - (100*(a + c - Cos[(-5*a + 5*c)^2]/100))/201


Now we just need to figure out how to apply this to subsequences of whatever list you want to pass in. Here's a tricky/fancy way:

Partition[{x1, x2, x3, x4, x5, x6, x7, x8, x9}, 3, 1, {2, -2}, {1, Null, 0}, newF]


Note how we used some padding to make sure that a 0 was inserted at the beginning and a 1 at the end.

Here is a way that's probably clearer because the padding is more explicit:

BlockMap[Apply[newF], Join[{0}, {x1, x2, x3, x4, x5, x6, x7, x8, x9}, {1}], 3, 1]


You could put this into a function which could work on either a list as its argument or a sequence. There are probably several other ways you could restructure the list or apply the function to sublists. I haven't looked closely at the rest of your question, but I presume you could do a similar thing with g.

You can define a vector-valued pure function that returns $$\{f_1(x_1,x_2,\ldots,x_m),f_2(x_1,x_2,\ldots,x_m), \cdots, f_m(x_1,x_2,\ldots,x_m)\}$$ when input $$\{x_1,x_2,\ldots,x_m\}$$ (similarly for gs) to simplify the first 40+ lines of your code into 4 lines:

ClearAll[fvec, gvec]

fvec = # - 100/201 (Prepend[0] @ Most[#] + Append[1] @ Rest[#] -
1/100 Cos[(5 # - 5 #2)^2 & @@@
Transpose[{Append[1] @ Rest[#], Prepend[0]@Most[#]}]]) &;

gvec = # - fvec @ # &;


Examples:

vars = Array[Subscript[x, #] &, 3];

fvec @ vars  // Column


gvec @ vars // Column


jacobian = D[gvec @ vars, {vars}];

MatrixForm @ jacobian


vars = Array[Subscript[x, #] &, 9];

fvec @ vars  // Column


vars = Array[Subscript[x, #] &, 30];

fvec @ vars  // Column
`

etc..