2
$\begingroup$

I am interested in solving a problem that takes the following form. I solve, using some numerical method such as RK4, an ODE of the form

$$ r(v)L(v;\{a_i\})+r^{\prime\prime}(v)=0\,, $$

where $L(v;\{a_i\})$ is a function of $v$ and some parameters $\{a_i\}_{i=0,\ldots,N-1}$. I know the initial conditions $r(1)=1$ and $r^\prime(1)=0$.

For each set of $\{a_i\}$, I compute $r(v)$, together with an integral of the form

$$ Z(\{a_i\})=\int_0^1{\rm d}v\,r(v)^2B(v;\{a_i\})\,, $$

where, once again, $B(v;\{a_i\})$ depends on $v$ and on the $\{a_i\}$. I would like to find the $\{a_i\}$ that maximise $Z(\{a_i\})$, restricted to configurations for which $r(v)>0$. Any idea how to do this? I tried using the in-built routines NMaximize together with NDSolve, but could not find a way to tell NMaximize to restrict to solutions $r(v)$ for which $r(v)$ is positive. Ideally, I want to take $N$ to be rather large, like $N=100$.

For example, take $L(v;\{a_i\})=[\partial_v \Phi(v;\{a_i\})]^2$ and $B(v;\{a_i\})=\Phi(v;\{a_i\})^2$ with

$$ \Phi(v;\{a_i\})=\sum_{i=0}^{N-1}a_iv^i,\quad \text{where}\quad a_i\in\mathbb{R}\,. $$

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ Could you show expressions for L and B? $\endgroup$ Commented Oct 1 at 15:33
  • 2
    $\begingroup$ Not nearly enough information to allow readers to provide assistance. $\endgroup$ Commented Oct 1 at 15:54
  • $\begingroup$ Thanks for the feedback. I just added an example. $\endgroup$ Commented Oct 1 at 20:25
  • $\begingroup$ @user12588 Do you mean this $\Phi(v;\{a_i\})=\sum_{n=0}^Na_iv^i$? $\endgroup$ Commented Oct 2 at 3:42
  • $\begingroup$ Yes - it has been amended. $\endgroup$ Commented Oct 2 at 9:19

1 Answer 1

Reset to default
2
$\begingroup$

This is toy example to demonstrate that this inverse scattering problem is unstably with the given initial data. First, we compute solution at $N=12$ as follows

nn = 12; A = Array[a, nn]; V = Table[(x)^i, {i, 0, nn - 1}]; V1 = 
 Table[i (x)^(i - 1), {i, 0, nn - 1}]; Phi = A.V; Phi1 = A.V1; L = 
 Phi1^2; B = Phi^2;

asol = AsymptoticDSolveValue[{r[x] (Phi1^2) + r''[x] == 0}, 
  r[x], {x, 0, nn}];
asol1 = D[asol, x] /. x -> 1;
Needs["NumericalDifferentialEquationAnalysis`"]

gqw = GaussianQuadratureWeights[nn, 0, 1]; points = 
 gqw[[All, 1]]; weights = gqw[[All, 2]];

int = weights.(asol^2 B /. x -> points) // Chop;

cons = asol /. x -> points // Chop;

sol = NMaximize[{int, 
   Join[Table[cons[[i]] > 0, {i, nn}], {asol == 1, asol1 == 0} /. 
     x -> 1]}, Join[A, {C[1], C[2]}]];

This is solution at $N=12$

{0.341698, {a[1] -> -0.854314, a[2] -> 0.454189, a[3] -> 0.103341, 
  a[4] -> -0.333199, a[5] -> -0.0923221, a[6] -> 0.0244073, 
  a[7] -> -0.495512, a[8] -> -0.266041, a[9] -> 0.355511, 
  a[10] -> 0.0456651, a[11] -> -0.305091, a[12] -> 0.420078, 
  C[1] -> 0.491838, C[2] -> 0.577452}}

Visualization

pl = Plot[asol /. sol[[2]], {x, 0, 1}, Frame -> True, 
  PlotRange -> {0, 1}, FrameLabel -> {"v", "r"}, 
  PlotStyle -> {Red, Dashed}]

Figure 1

Second, we compute solution at $N=13$ using code above at nn=13, we have new solution

{1955.15, {a[1] -> -6.84304, a[2] -> -0.304688, a[3] -> -1.66847, 
  a[4] -> -0.0229783, a[5] -> 0.543985, a[6] -> 0.30122, 
  a[7] -> -0.648714, a[8] -> -0.522593, a[9] -> -1.25001, 
  a[10] -> -0.537674, a[11] -> -0.246614, a[12] -> 0.599196, 
  a[13] -> 0.665486, C[1] -> 6.5781, C[2] -> 2.11588}}

Visualization

Figure 2

What can we recommend in this regard?

  1. Constraints on the coefficients $a_i$ must be introduced to ensure the convergence of the series $\Phi$.
  2. BVP is more logic in this case.

In any case, this example demonstrates how to solve this problem under the given constraints $r(v)>0$.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.