1
$\begingroup$

I have the following ODE from the paper https://sci-hub.hkvisa.net/10.1016/j.aml.2018.02.016

$y^{\prime\prime}(t)=-ay^{\prime m}(t)+y(t)(b y^{2}(t)-\frac{3}{2}y(t)+\frac{1}{2})$. The boundary conditions are $y(0)=1$ and $y(1)=2$. When $m=2$ and $b=a=1$. The author obtained the attached figure by optimizing the iteration enter image description here

The details of this question is related with How to optimize the $L^{2}$ norm in Mathematica for an iteration scheme and How to plot its 3D graph?

$\endgroup$
1
  • 1
    $\begingroup$ 1. Include the equations as Mathematica code as well, particularly if this is a continuation of another question. 2) include your question explicitly and clearly in the text of the post; don't leave it only in the title. $\endgroup$
    – MarcoB
    Apr 3 at 14:08

1 Answer 1

3
$\begingroup$

We can use code from my answer here to plot the $L^2$-norm of the residual error and compute optimized parameters $\alpha, \beta$. Please, note, that in my code $a=\alpha, b= \beta$, don't mess up these parameters with a=1,b=1 in ODE.

nmax = 10; var = Table[t^n, {n, 0, nmax}];
nn = 1; x[0] = 1 + t;
Do[xs[n] = x[n] /. t -> s;
   gn = D[xs[n], {s, 2}] + D[xs[n], s]^2 - xs[n]^3 + 3/2 xs[n]^2 - 
     1/2 xs[n]; lst = CoefficientList[gn, t]; 
   g[n] = If[Length[lst] < Length[var], gn, 
     Take[lst, Length[var]] . var];
   in1 = Integrate[s*g[n], s]; in2 = Integrate[g[n], s] - in1;
   int1 = in1 /. s -> t; int2 = (in2 /. s -> 1) - (in2 /. s -> t);
   yn = x[n] + b*(1 - t) int1 + b*t int2; ys[n] = yn /. t -> s; 
   fn = D[ys[n], {s, 2}] + D[ys[n], s]^2 - ys[n]^3 + 3/2 ys[n]^2 - 
     1/2 ys[n]; lst = CoefficientList[gn, t]; 
   f[n] = If[Length[lst] < Length[var], fn, 
     Take[lst, Length[var]] . var]; inf1 = Integrate[s*f[n], s]; 
   inf2 = Integrate[f[n], s] - inf1;
   intf1 = inf1 /. s -> t; intf2 = (inf2 /. s -> 1) - (inf2 /. s -> t);
   xn = x[n] + a (1 - t) intf1 + a*t*intf2;
   lst = CoefficientList[xn, t];
   x[n + 1] = 
    If[Length[lst] < Length[var], xn, 
     Take[lst, Length[var]] . var];, {n, 0, nn}] // Quiet;

We can't compute the $L^2$-norm of the residual error directly for this case, and instead of it we use numerical algorithm

cl = CoefficientList[(1 - t) intf1 + t*intf2, t] // Simplify;
res = Take[cl, Length[var]] . var // N;
error1 = 
  Table[{a1, b1, 
    a1^2 NIntegrate[res^2 /. {a -> a1, b -> b1}, {t, 0, 1}]}, {a1, 1, 
    1.6, .1}, {b1, .0, .6, .1}];

 er1 = Interpolation[Flatten[error1, 1]];

We use function er1 to minimize error

NMinimize[{er1[a, b], 1 <= a <= 1.6, .0 <= b <= .6}, {a, b}]

(*Out[]= {2.50304*10^-6, {a -> 1.45582, b -> 0.378981}}*)

Visualization

Plot3D[er1[a, b], {a, 1.2, 1.55}, {b, 0.3, .5}, 
 ColorFunction -> Hue, AxesLabel -> {"\[Alpha]", "\[Beta]", ""}, 
 Boxed -> False, PlotTheme -> "Marketing"]

Figure 1

Since parameters $\alpha=1.45582, \beta= 0.378981$ are very differ from that in the paper, α = 0.9345414427 and β = 0.8817524743, and also we don't know exact solution we need to verifier numerical solution with other method. First we compute x[n+1] with our code as follows

nmax = 60; var = Table[t^n, {n, 0, nmax}];
nn = 3; x[0] = 1 + t; a = 1.4558248634975057; b = 0.3789811601409471;
Do[xs[n] = x[n] /. t -> s;
   gn = D[xs[n], {s, 2}] + D[xs[n], s]^2 - xs[n]^3 + 3/2 xs[n]^2 - 
     1/2 xs[n]; lst = CoefficientList[gn, t]; 
   g[n] = If[Length[lst] < Length[var], gn, 
     Take[lst, Length[var]] . var];
   in1 = Integrate[s*g[n], s]; in2 = Integrate[g[n], s] - in1;
   int1 = in1 /. s -> t; int2 = (in2 /. s -> 1) - (in2 /. s -> t);
   yn = x[n] + b*(1 - t) int1 + b*t int2; ys[n] = yn /. t -> s; 
   fn = D[ys[n], {s, 2}] + D[ys[n], s]^2 - ys[n]^3 + 3/2 ys[n]^2 - 
     1/2 ys[n]; lst = CoefficientList[gn, t]; 
   f[n] = If[Length[lst] < Length[var], fn, 
     Take[lst, Length[var]] . var]; inf1 = Integrate[s*f[n], s]; 
   inf2 = Integrate[f[n], s] - inf1;
   intf1 = inf1 /. s -> t; intf2 = (inf2 /. s -> 1) - (inf2 /. s -> t);
   xn = x[n] + a (1 - t) intf1 + a*t*intf2;
   lst = CoefficientList[xn, t];
   x[n + 1] = 
    If[Length[lst] < Length[var], xn, 
     Take[lst, Length[var]] . var];, {n, 0, nn}] // Quiet; 

Then we prepare list for comparision

y= Table[{t, x[3]}, {t, 0, 1, .01}]

(*Out[]= {{0., 1.}, {0.01, 1.01241}, {0.02, 1.02466}, {0.03, 1.03677}, {0.04,

 
  1.04874}, {0.05, 1.06057}, {0.06, 1.07226}, {0.07, 1.08382}, {0.08, 
  1.09526}, {0.09, 1.10657}, {0.1, 1.11777}, {0.11, 1.12885}, {0.12, 
  1.13981}, {0.13, 1.15067}, {0.14, 1.16142}, {0.15, 1.17207}, {0.16, 
  1.18262}, {0.17, 1.19308}, {0.18, 1.20344}, {0.19, 1.21371}, {0.2, 
  1.2239}, {0.21, 1.23401}, {0.22, 1.24403}, {0.23, 1.25398}, {0.24, 
  1.26385}, {0.25, 1.27365}, {0.26, 1.28338}, {0.27, 1.29305}, {0.28, 
  1.30265}, {0.29, 1.3122}, {0.3, 1.32168}, {0.31, 1.33111}, {0.32, 
  1.34049}, {0.33, 1.34983}, {0.34, 1.35911}, {0.35, 1.36835}, {0.36, 
  1.37755}, {0.37, 1.38671}, {0.38, 1.39583}, {0.39, 1.40493}, {0.4, 
  1.41399}, {0.41, 1.42302}, {0.42, 1.43202}, {0.43, 1.44101}, {0.44, 
  1.44997}, {0.45, 1.45891}, {0.46, 1.46784}, {0.47, 1.47675}, {0.48, 
  1.48566}, {0.49, 1.49456}, {0.5, 1.50345}, {0.51, 1.51233}, {0.52, 
  1.52122}, {0.53, 1.53011}, {0.54, 1.539}, {0.55, 1.54791}, {0.56, 
  1.55682}, {0.57, 1.56574}, {0.58, 1.57468}, {0.59, 1.58363}, {0.6, 
  1.59261}, {0.61, 1.6016}, {0.62, 1.61062}, {0.63, 1.61967}, {0.64, 
  1.62875}, {0.65, 1.63787}, {0.66, 1.64701}, {0.67, 1.6562}, {0.68, 
  1.66543}, {0.69, 1.6747}, {0.7, 1.68402}, {0.71, 1.69338}, {0.72, 
  1.7028}, {0.73, 1.71228}, {0.74, 1.72181}, {0.75, 1.7314}, {0.76, 
  1.74105}, {0.77, 1.75077}, {0.78, 1.76057}, {0.79, 1.77043}, {0.8, 
  1.78037}, {0.81, 1.79038}, {0.82, 1.80048}, {0.83, 1.81066}, {0.84, 
  1.82093}, {0.85, 1.83129}, {0.86, 1.84175}, {0.87, 1.8523}, {0.88, 
  1.86296}, {0.89, 1.87372}, {0.9, 1.88458}, {0.91, 1.89556}, {0.92, 
  1.90665}, {0.93, 1.91787}, {0.94, 1.9292}, {0.95, 1.94066}, {0.96, 
  1.95225}, {0.97, 1.96398}, {0.98, 1.97584}, {0.99, 1.98785}, {1., 
  2.}}*)

Third we compute numerical solution with using collocation method and Euler wavelets

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 = 3; M0 = 4; nn = 
 Total[With[{k = k0, M = M0}, 
   Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; 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]; Psi[y_] := Psijk /. t1 -> y; 
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;
M = nn; var = Array[v, {M}];
X2[t_] := var . Psi[t]; X1[t_] := var . int1[t] + v1; 
X[t_] := var . int2[t] + v1 t + v0;
eqn = Join[
  Table[X2[s] + X1[s]^2 - X[s]^3 + 3/2 X[s]^2 - 1/2 X[s] == 0, {s, 
    tcol}], {X[0] == 1, X[1] == 2}]; varM = Join[{v0, v1}, var];

sol = FindRoot[eqn, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]  

Finally we plot two solutions and have a good coincidence

Show[Plot[X[t] /. sol, {t, 0, 1}, AxesLabel -> {"t", "y"}], 
 ListPlot[y, PlotStyle -> Red]] 

Figure 2

Third numerical solution we can compute with

ns = NDSolve[{z''[s] + z'[s]^2 - z[s]^3 + 3/2 z[s]^2 - 1/2 z[s] == 0, 
   z[0] == 1, z[1] == 2}, z, {s, 0, 1}]

Differences between 3 solutions shown below.

{ListPlot[
  Table[X[s] /. sol, {s, 0, 1, 0.01}] - 
   Table[z[s] /. ns[[1]], {s, 0, 1, 0.01}], PlotRange -> All, 
  PlotLabel -> "Collocation - NDSolve"], 
 ListPlot[Table[X[s] /. sol, {s, 0, 1, 0.01}] - y[[All, 2]], 
  PlotRange -> All, PlotLabel -> "Collocation - PredictorCorrector"], 
 ListPlot[Table[z[s] /. ns[[1]], {s, 0, 1, 0.01}] - y[[All, 2]], 
  PlotRange -> All, PlotLabel -> "NDSolve - PredictorCorrector"]}

Figure 3

We can decrease error of colocation method up to $10^{-13}$ then the difference between NDSolve solution and colocation method is about $2\times 10^{-8}$ - see Figure 4 below. We also can compute predictor-corrector solution with α = 0.9345414427 and β = 0.8817524743 from the paper. In this case the error increases from $10^{-5}$ (see Figure 3 above) to $4\times 10^{-2}$ - see Figure 4 below. It means that parameters of the predictor-corrector algorithm have not been optimized in the paper cited.
Figure 4

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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