# Taylor approximation of integrals

I have BVP as in the picture and I want to solve it numerically by method given as follow. In the original paper, it is motioned that first iteration is computed and other are computed with Taylor approximation of the integral inside the BVP. But I do not How I replace this integral by Taylor or any other so Mathematica Show numerical values fast.

\[Delta] = 10^-20;
Clear[x];
x[0 _] := x[0] = Function[t, 1.000000];
a[n_] := a[n] = 0.99;
x[n_] := x[n] = Function[t, Evaluate[Chop[Expand[x[n - 1][t] + Integrate[Expand[((t^3/3 - t^2/2 + 1/6) s^3 - (t^3/2 - t^2 + t/2) s^2) (x[n - 1]''''[s] - (0.65 Integrate[(x[n - 1]'[s])^2, {s, 0, 1}]) + 0.20/(x[n - 1][s])^4 + 5/(-0.396 + x[n - 1][s])^2 + 0.325/x[n - 1][s])], {s, 0, t}] + Integrate[Expand[((t^3/3 - t^2/2) s^3 - (t^3/2 - t^2) s^2 - (t^2*s)/2 + t^3/6) (x[n - 1]''''[s] - (0.65 Integrate[(x[n - 1]'[s])^2, {s, 0, 1}]) + 0.20/(x[n - 1][s])^4 + 5/(-0.396 + x[n - 1][s])^2 + 0.325/x[n - 1][s])], {s, t, 1}]], \[Delta]]]]
a1a = Table[x[n][0.1], {n, 0, 1}]


Where $$x_{0}(t)=1.000000$$ is the initial iterate. When I increase the number of iterations from 1 to 6 or higher (i.e. {n,0,6}) then Mathematica doest not show results fast. Now I want to replace the integral inside the code i.e.

Integrate[(x[n - 1]'[s])^2, {s, 0, 1}]


by Taylor approximation or any other as suggested in the paper. full paper is here https://doi.org/10.1007/s10958-022-05886-w

• Do you mean equation (2.1) with boundary conditions (2.2) in the paper A NUMERICAL INVESTIGATION OF THE BUCKLING OF DOUBLY CLAMPED NANO-ACTUATORS GOVERNED BY AN INTEGRO-DIFFERENTIAL EQUATION? Dec 1, 2023 at 4:33

I found one scheme that converges, but with not very high accuracy. Note that increasing the number of terms of the series and the number of iterations does not reduce the error. It is possible that this code can be used for parametric studies, although it is not very fast.

nmax = 19; niter = 20; alpha = 20; beta = 5; gamma = 0.325; eta = 0.96; k = -0.396; P = 0;
Clear[x];
x[0] = 1; int01 = 0;
Do[x0 = x[n]; x4 = D[x0, {t, 4}]; x2 = D[x0, {t, 2}];
int0 = Normal[
Series[x4 - (eta int01 + P) x2 + alpha/(x0)^4 + beta/(k + x0)^2 +
gamma/x0, {t, 1/2, nmax}]] /. t -> s;
x[n + 1] =
x[n] + Integrate[((t^3/3 - t^2/2 + 1/6) s^3 - (t^3/2 - t^2 +
t/2) s^2) (int0), {s, 0, t}] +
Integrate[((t^3/3 - t^2/2) s^3 - (t^3/2 - t^2) s^2 - (t^2*s)/2 +
t^3/6) (int0), {s, t, 1}]; x1 = D[x[n + 1], {t, 1}];
int01 = NIntegrate[(x1)^2, {t, 0, 1}];, {n, 0,
niter}]


It takes about 60s, and finally we have for maximum deflection

Table[{n, 1 - x[n] /. t -> .5}, {n, 1, niter + 1}]

Out[]= {{1, 0.0886212}, {2, 0.116071}, {3, 0.126898}, {4,
0.13152}, {5, 0.133545}, {6, 0.134441}, {7, 0.134838}, {8,
0.135015}, {9, 0.135093}, {10, 0.135128}, {11, 0.135143}, {12,
0.13515}, {13, 0.135153}, {14, 0.135155}, {15, 0.135155}, {16,
0.135155}, {17, 0.135156}, {18, 0.135156}, {19, 0.135156}, {20,
0.135156}, {21, 0.135156}}


We can compare this result with Table 3 from the paper (note, that legend to this table has a typo with parameters)

Our proposed PGS 0.135526676798
Numerical [30] 0.13536


We can check their more accurate result 0.135526676798 using our algorithm from the paper. First line

Get["NumericalDifferentialEquationAnalysis"]


Then

UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2)  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 = 4; M0 = 7; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
tcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; 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]; Int3 = Integrate[Int2, t1]; Int4 =
Integrate[Int3, t1]; Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;
int3[y_] := Int3 /. t1 -> y; int4[y_] := Int4 /. t1 -> y;
X = Array[xx, {nn}]; u4[y_] := X . Psi[y]; u3[y_] := X . int1[y] + p0;
u2[y_] := X . int2[y] + p0 y + p1;
u1[y_] := X . int3[y] + p0 y^2/2 + p1 y + p2;
u0[y_] := X . int4[y] + p0 y^3/6 + p1 y^2/2 + p2 y + p3;
np = 100; g = GaussianQuadratureWeights[np, 0, 1];
gaussint[f_, z_] := Sum[g[[i, 2]] f /. z -> g[[i, 1]], {i, np}]; int =
Sum[g[[i, 2]] u1[z]^2 /. z -> g[[i, 1]], {i, np}];

With[{eta = 0.96, P = 0, alpha = 20, beta = 5, gamma = 0.325,
kappa = -0.396},
eq = (u4[s] - (eta int + P)*u2[s] + alpha/(u0[s])^4 +
beta/(kappa + u0[s])^2 +
gamma/u0[s])]; With[{eps = \$MachineEpsilon},
bc = {u0[eps] == 1, u0[1 - eps] == 1, u1[eps] == 0,
u1[1 - eps] == 0}];
eqn = Table[eq == 0, {s, xcol}];
var = Join[X, {p0, p1, p2, p3}];
sol = FindRoot[Join[eqn, bc],
Table[{var[[i]], RandomReal[]}, {i, Length[var]}],
MaxIterations -> 1000];


Finally we have for maximum deflection

1 - u0[1/2] /. sol

Out[]= 0.13552748281600993


Therefore our result agrees with 0.135526676798 with error of 10^-6. But we don't know how they compute this result. We can improve our computation by doubling points np=200,k0 = 5; M0 = 7; so nn=112, and maximum deflection is 0.13552748274091223. We also can plot function u[x] as

Plot[u0[s] /. sol, {s, 0, 1}]


Note, that we reproduce Figure 1 from the paper.

There are some typos in your equations if what you want to reproduce is table 2 and figure 1. In particular the alpha parameter is 20, not 0.2, the nu parameter is .96, not .65, and there is a second derivative to be multiplied bu those inner integrals.

For speed I made several other changes.

(1) Each x[n+1] is computed from x[n] as an interpolating function. I used fairly high order of 11 because we need continuity and there are fourth derivatives in the integrands.

(2) I use numerical integration.

(3) The inner integrals are constant for each n so, for given n, I compute them one time and use those values in (what were) the outer integrals.

(4) I was not able to see where the value xi was actually used. The authros have it set to 4, but...

When I break the unit interval into hundredths each iteration after x[1] takes around 17 seconds on my desktop machine. Not sure if the results would converge faster with finer granularity.

ClearAll[x];
x[0] = Function[t, 1];
alpha = 20;
beta = 5;
gamma = .325;
kappa = -.396;
nu = .96;
xi = 4;
x[n_] := x[n] = Module[{vals, int},
vals = Table[
int = nu* NIntegrate[(x[n - 1]'[s])^2, {s, 0, 1}];
{{t},
x[n - 1][t] +
NIntegrate[((t^3/3 - t^2/2 + 1/6) s^3 - (t^3/2 - t^2 +
t/2) s^2) (x[n - 1]''''[s] - int*x[n - 1]''[s] +
alpha/(x[n - 1][s])^4 + beta/(kappa + x[n - 1][s])^2 +
gamma/x[n - 1][s]), {s, 0, t}] +
NIntegrate[((t^3/3 - t^2/2) s^3 - (t^3/2 - t^2) s^2 - (t^2*s)/
2 + t^3/6) (x[n - 1]''''[s] - int*x[n - 1]''[s] +
alpha/(x[n - 1][s])^4 + beta/(kappa + x[n - 1][s])^2 +
gamma/x[n - 1][s]), {s, t, 1}]}, {t, 0., 1, .01}];
Interpolation[vals, InterpolationOrder -> 11]]


Here are some iterations.

Timing[plots = Table[Quiet[x[n]];
Plot[x[n][t], {t, .01, .99}], {n, 1, 10}]]


They quickly become indistinguishable to the eye. We can assess convergence by plotting differences.

Table[Plot[x[n][t] - x[n - 1][t], {t, .01, .99}], {n, 1, 10}]


They seem to drop by about a factor of 2 with each iteration. Which is in rough agreement with table 2 and related remarks in the paper.

--- edit ---

If I change the interpolation point increment from .01 to .05 this become much faster, taking under 8 seconds for 10 iterations. And the convergence is about the same as before. After 20 iterations the maximum difference between the last pair is about 10^(-8).

--- end edit ---

• This is a good approach (+1). But how can we be sure that this is correct? The quoted article has several typos. For example, the data in Table 3 clearly does not correspond to what is stated there. If you use your algorithm and compute with their parameters, alpha = 25; beta = 15; gamma = 0.65; kappa =0; nu = 0; xi = 4; you will not get their answer. They also refer to paper [30], but in [30] were used different parameters alpha = 20; beta = 5; gamma =?; kappa =?; nu = 0; xi = 4; shown in Table 1. Dec 2, 2023 at 5:21
• @AlexTrounev I’ll take a look when I’m back at my computer. But offhand the issue could be at their end (e.g mismatches between claimed and actual parameters used) or, more likely, an issue with the numerics in my code. Dec 2, 2023 at 16:46
• @AlexTrounev Some trial-and-error shows that their Table 3 has the correct values for the claimed parameter settings EXCEPT with beta=0 rather than the stated value of 15. Dec 4, 2023 at 20:23
• Please, see my answer where used set {eta = 0.96, P = 0, alpha = 20, beta = 5, gamma = 0.325, kappa = -0.396} to compute maximum deflection 0.135526676798 shown in Table 3. Note, that in my answer this value computed with 2 different codes, one is some implementation of method proposed in the paper cited, and another is the Euler wavelet collocation method. It is why I think that legend to the Table 3 has a typo. Unfortunately, I was not able to calculate this value using your code. If you can do it yourself, please show me how. Dec 5, 2023 at 8:23
• @AlexTrounev I fixed a problem in my code wherein I had copied the version using the Mann rather than Picard scheme for iterations. The former is just not stable for my method. With that change, and with setting nu=1 and other parameters as you provide in your comment, I get convergence of the deflection at .5 to the value 0.135504. Dec 5, 2023 at 19:27