9
$\begingroup$

I want to solve a second order differential equation in the interval[-1:1], which does not have a analytic solution, \begin{eqnarray} y''(x) &=& k \phi^2(x)y(x) \\ \phi(x) &=& \frac{1}{2}\left(1-\sin\left[ \frac{\pi}{2}x\right] \right)\\ y'[1] &=& 1 \\ y[1] &=& 1 \end{eqnarray}

It is possible to get numerical solution of the equation using NDSolve. But I want to approach the solution iteratively starting with a Initial guess function for example, $y_0(x) = \frac{1}{2}(1+x)$ with appropriate boundary conditions, \begin{eqnarray} y_{n+1}''(x) &=& k\phi^2(x)y_n(x) \end{eqnarray}

To do this I have following code,

ClearAll;
k = 3955/100;
W = 1
phi[x_] := (1/2)*(1 - Sin[(Pi/2)*x/W])
vE[x_] := Piecewise[{{0, -W <= x < 0}, {x/W, 0 <= x <= W}}];
pcw[x_] := (W + x)/(2*W);
so = NDSolve[{u''[x] == k*phi[x]*phi[x]*u[x], u[-W] == 0, u[W] == 1}, 
  u, {x, -W, W}, WorkingPrecision -> 22, InterpolationOrder -> All]
dd[x_] = Q[x] /. 
   First@DSolve[{Q''[x] == k*phi[x]*phi[x]*pcw[x], Q[W] == 1, 
      Q'[W] == 1}, Q, x];
Plot[Evaluate[{dd[t], u[t] /. so, vE[t], pcw[t]}], {t, -W, W}]
dd[x_] = Q[x] /. 
   First@DSolve[{Q''[x] == k*phi[x]*phi[x]*dd[x], Q[W] == 1, 
      Q'[W] == 1}, Q, x];
Plot[Evaluate[{dd[t], u[t] /. so, vE[t], pcw[t]}], {t, -W, W}]
dd[x_] = Q[x] /. 
   First@DSolve[{Q''[x] == k*phi[x]*phi[x]*dd[x], Q[W] == 1, 
      Q'[W] == 1}, Q, x];
Plot[Evaluate[{dd[t], u[t] /. so, vE[t], pcw[t]}], {t, -W, W}]
dd[x_] = Q[x] /. 
   First@DSolve[{Q''[x] == k*phi[x]*phi[x]*dd[x], Q[W] == 1, 
      Q'[W] == 1}, Q, x];
Plot[Evaluate[{dd[t], u[t] /. so, vE[t], pcw[t]}], {t, -W, W}]

Which is certainly not elegant way getting things done and also it is taking too long to run than I was expecting. How can I use recursion in this case?

The numeric solution is:

Plot[u[x] /. so, {x, -W, W}]

Mathematica graphics

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ 1) ClearAll; accomplishes nothing; maybe you want Clear["Global*"]. 2) Each of your definitions of pcw` overrides the previous one, so that only the last definition pcw[x_] := (W + x)/(2*W) remains in effect at the end. Did you intend that? $\endgroup$
    – m_goldberg
    Commented Oct 26, 2014 at 7:07
  • $\begingroup$ Thank you m_goldberg for the comments..! I will use ClearAll properly now. for the pcw function any definition is good enough. $\endgroup$
    – chatur
    Commented Oct 27, 2014 at 9:21
  • 2
    $\begingroup$ If you simply want to make the code elegant, have a look at Do, While, Nest, NestWhile etc. Nevertheless, I'm afraid it's hard to speed up the code if you want dd[x] to be analytic. BTW, shouldn't the phi[x] phi[t] be phi[x]^2? $\endgroup$
    – xzczd
    Commented Oct 28, 2014 at 2:53
  • $\begingroup$ @xzczd, thanks for pointing out that. phi[x]*phi[t] should be phi[x]^2 $\endgroup$
    – chatur
    Commented Oct 28, 2014 at 9:36
  • $\begingroup$ I guess you want an iteration method as in the proof of e.g. Picard's theorem rather thant recurssion. Unfortunately I have no time for answering any questions. $\endgroup$
    – Artes
    Commented Oct 28, 2014 at 9:41

1 Answer 1

Reset to default
3
$\begingroup$

The following is perhaps elegant but not that fast. The main problem is that your ODE is increasing in complexity at each iteration. Applying N[...] at each step helps with that but sooner or later (fatally) you'll lose precision and the thing will stop converging to the desired solution.

k = 3955/100;
w = 1;
phi[x_] := (1/2) (1 - Sin[Pi/2 x/w])
pcw[x_] := (w + x)/(2 w);
t[f_]:= (q@x/.First[DSolve[{q''@x== k phi@x^2 f, q@w== 1, q'@w== 1}, q, x]]//Simplify) 

nl = NestList[ t[#] &, t[pcw@x], 5];
(* The numeric solution*)
nds = NDSolveValue[{y''[x] == k phi[x]^2 y[x], y'[1] == 1, y[1] == 1},y, {x, -1, 1}]

Show[Plot[nds[x], {x, -1, 1}, PlotStyle -> {Thick, Red}], 
     Plot[nl, {x, -w, w}], PlotRange -> All]

Mathematica graphics

here you can see the increasing complexity

LeafCount /@ nl
(* {88, 312, 797, 1659, 3011, 4965} *)
Improve this answer
$\endgroup$

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