# FindMinimum running forever

I cannot get any output from the below-mentioned optimization problem. The program keeps running forever. What can I do to solve this optimization problem successfully?

ClearAll["Global*"];
b=19;
a=15;
es=1;                                                 (*Element size for FEA*)
noofnode=IntegerPart[((b-a)/es)+1];
init=0&;
eb=380*10^9;     (*MPa*)
ea=71*10^9;       (*MPa*)
new=0.3;                             (*Poisson's Ratio*)
alphab=8.0*10^(-6);
alphaa=23.1*10^(-6);
rowb=0.96;
rowa=2.70;
ta=20;                      (*Degree C*)
tb=150;
rpm=150.0;
wom=(2.0*\[Pi]*rpm)/60;
kb=35.0;
ka=204.0;
lam=(1/(b-a))*Log[kb/ka];               (*Constant for Exponential Temp Profile*)
f[ri_]:=ExpIntegralEi[-lam*ri];
c1=(tb-ta)/(f[b]-f[a]);
c2=((ta*f[b])-(tb*f[a]))/(f[b]-f[a]);
T[r_]:=c2+c1*ExpIntegralEi[-lam*r];
ed[v_]:=(v*ea)+((1-v)*eb);
rowd[v_]:=(v*rowa)+((1-v)*rowb);
e[v1_,v2_,v3_,v4_,v5_]:=Interpolation[{ed[v1],ed[v2],ed[v3],ed[v4],ed[v5]}];
row[v1_,v2_,v3_,v4_,v5_]:=Interpolation[{rowd[v1],rowd[v2],rowd[v3],rowd[v4],rowd[v5]}];
z[v1_?NumericQ,v2_?NumericQ,v3_?NumericQ,v4_?NumericQ,v5_?NumericQ]:=NDSolveValue[{(3+new) r^2 wom^2 row[v1,v2,v3,v4,v5][r-14]+F'[r]+(e[v1,v2,v3,v4,v5]'[r-14] (new F[r]-r (r^2 wom^2 row[v1,v2,v3,v4,v5][r-14]+F'[r])))/e[v1,v2,v3,v4,v5][r-14]+r^3 wom^2 row[v1,v2,v3,v4,v5]'[r-14]+r e[v1,v2,v3,v4,v5][r-14] (T[r] alpha[v1,v2,v3,v4,v5]'[r-14]+alpha[v1,v2,v3,v4,v5][r-14] T'[r])+r F''[r]==F[r]/r,F[a]==0,F[b]==0},F,{r,a,b}];
o[v1_,v2_,v3_,v4_,v5_]:=Sqrt[Sum[(((z[v1,v2,v3,v4,v5][j])*10^(-9))^2),{j,15,19}]/5];
FindMinimum[{o[v1,v2,v3,v4,v5],{0<=v1<=1,0<=v2<=1,0<=v3<=1,0<=v4<=1,0<=v5<=1}},{v1,v2,v3,v4,v5}]

• The problem is that NDSolve can't find solution for some combination of parameters. Also FindMinimum with NDSolve is very slow combination. It could be better to optimize some FDM solution. Commented Feb 5, 2022 at 13:35
• I don't understand what equation do you try to solve with this method? Can you show equation and boundary conditions? Commented Feb 6, 2022 at 9:44
• @AlexTrounev I have edited the post by adding your requirements. Please go through the edited post. Commented Feb 6, 2022 at 12:56
• Is the number of nodal points is fixed (5 in your code) or it could be 6, 7, 8...? Commented Feb 6, 2022 at 13:31

This problem can be solved with using the Euler wavelets collocation method explained in our paper and on the page. First we map region $$a\le r\le b$$ on the unit interval by substitution r=a+(b-a) x in GE. We also define F2=F'',F1=F',v1=v', then we have

ClearAll["Global*"];
b = 19;
a = 15;
eb = 380*10^9;(*MPa*)ea =
71*10^9;(*MPa*)new = 0.3;(*Poisson's Ratio*)alphab = 8.0*10^(-6);
alphaa = 23.1*10^(-6);
rowb = 0.96;
rowa = 2.70;
ta = 20;(*Degree C*)tb = 150;
rpm = 150.0;
wom = (2.0*\[Pi]*rpm)/60;
kb = 35.0;
ka = 204.0;
L = b - a;
lam = (1/(b - a))*Log[kb/ka];(*Constant for Exponential Temp Profile*)
f[ri_] := ExpIntegralEi[-lam*ri];
c1 = (tb - ta)/(f[b] - f[a]);
c2 = ((ta*f[b]) - (tb*f[a]))/(f[b] - f[a]);
T = c2 + c1*ExpIntegralEi[-lam*r]; T1 = c1 E^(-lam r)/r;
ed = (v*ea) + ((1 - v)*eb);
alphad = (v*alphaa) + ((1 - v)*alphab);
rowd = (v*rowa) + ((1 - v)*rowb);
e1 = D[ed, v] v1[x]/L;
row1 = D[rowd, v] v1[x]/L;

eq = {(3 + new) r^2 wom^2 rowd +
F1[x]/L + (e1 (new F[x] - r (r^2 wom^2 rowd + F1[x]/L)))/ed +
r^3 wom^2 row1 + r ed (T alpha1 + alphad T1/L) + r F2[x]/L^2 -
F[x]/r, F[a] == 0, F[b] == 0} /. {v -> v[x], r -> a + L x,
F[a] -> F[0], F[b] -> F[1]} // Simplify;


Second, we define wavelets and all functions to be computed

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 3; With[{k = k0, M = M0},
var0 = Flatten[Table[cvar[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]];
nn = Length[var0];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; zcol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2,
nn + 1}]; ycol = zcol; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1]; Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;
M = nn; U1 = Array[u1, {M}]; U2 = Array[u2, {M}];
F2[x_] := U1 . Psi[x]; F1[x_] := U1 . int1[x] + u0;
F[x_] := U1 . int2[x] + u0 x + u00; v1[x_] := U2 . Psi[x];
v[x_] := U2 . int1[x] + v0;


Now we can solve GE and define F as function of v in collocation points

var = Join[U1, {u0, u00}]; eq1 =
Join[Table[eq[[1]] == 0, {x, xcol}], {eq[[2]], eq[[3]]}];

sol=Solve[eq1, var];

opt = Table[F[x] /. sol[[1]], {x, xcol}];


Finally we solve optimization problem with constraints $$0\le v\le 1$$ as follows

var2 = Join[U2, {v0}]; res =
Table[{v[x] >= 0, v[x] <= 1}, {x, Join[{0}, xcol, {1}]}]//Flatten; sol1 = NMinimize[{opt . opt, res}, var2];


Visualization

    {Plot[Evaluate[v[(r - a)/L] /. sol1[[2]]], {r, a, b},
PlotRange -> All, Frame -> True, PlotLabel -> "v"],
Plot[Evaluate[F[(r - a)/L] /. sol[[1]] /. sol1[[2]]], {r, a, b},
PlotRange -> All, Frame -> True, PlotLabel -> "F"]}


Update 1. We can force computation with option Method ->"DifferentialEvolution", for instance, for k0 = 2; M0 = 4;we define

scale =
Evaluate[opt[[1]] /.
Table[var2[[i]] -> RandomReal[{-1, 1}], {i, Length[var2]}]]

(*Out[]= -2.28631*10^7*)

sol1 =
NMinimize[{Norm[opt/scale], res}, var2,
Method -> "DifferentialEvolution"]

(*Out[]= {0.0000203025, {u2[1] -> -0.399257, u2[2] -> 0.107406,
u2[3] -> -0.0285942, u2[4] -> 0.0115987, u2[5] -> -0.267057,
u2[6] -> 0.042192, u2[7] -> -0.00704157, u2[8] -> 0.00169105,
v0 -> 1.}}*)


Visualization

• e1 = D[e, v] v1[x]/L; alpha1 = D[alphad, v] v1[x]/L; row1 = D[rowd, v] v1[x]/L; Can you please explain these lines? I think that should be ed instead of e. What is v1[x] for? and why did you divide by L in some places of the GE? In this model do e, alpha, and row are changing with the change of v. They must change with the change of v as I mentioned in the image. Finally, how much it took to give the mentioned outputs? Commented Feb 7, 2022 at 10:30
• It took less than a minute when e1 = D[e, v] v1[x]/L where e1 is giving 0. It should be e1 = D[ed, v] v1[x]/L. After putting e1 = D[ed, v] v1[x]/L the code is still running. Commented Feb 7, 2022 at 10:41
• @SazedurRahman Yes, you are right, it should be ed, Code has been updated. v[x], v1[x]=v'[x] used instead of interpolating in your code. e, alpha, row are changing with the change of v[x], and consequently e1, alpha1, row1 with v1[x]=v'[x]. Please, use updated version. Commented Feb 7, 2022 at 12:38
• I have another doubt. eq = {(3 + new) r^2 wom^2 rowd + F1[x]/L + (e1 (new F[x] - r (r^2 wom^2 rowd + F1[x]/L)))/ed + r^3 wom^2 row1 + r ed (T alpha1 + alphad T1/L) + r F2[x]/L^2 - F[x]/r, F[a] == 0, F[b] == 0} /. {v -> v[x], r -> a + L x, F[a] -> F[0], F[b] -> F[1]} // Simplify; Here why F1 and T are divided by L? Also, in e1 row1 alpha1 , v1 has been divided by L. Why all of these have been divided by L Commented Feb 7, 2022 at 13:57
• @SazedurRahman Please, note, that dr=L dx it is why all derivatives have 1/L. Commented Feb 7, 2022 at 14:14