Approximating ODE solution with polynomial function

Question

I am trying to decipher the method for generating non-linear frequency-modulated signals described in this paper by A. W. Doerry. On page 5, there is an algorithm described that boils down to approximating a solution to the following differential equation $$\omega'(t)=\frac{\omega'(0)}{W(\omega(t)-\omega(0))}$$ with the constraint of $$\int_{-T/2}^{+T/2}\omega'(t)\,dt=\Omega.$$ A few easy to work with forms for $\omega(t)$ are given. The one I'm focusing on is the polynomial form. $$\tilde\omega(t)=\sum_{n=1}^N c_n\, t^{n-1} $$ The suggested iterative procedure for finding the polynomial (or whichever form) function fit is

1. Start with some initial $\tilde\omega'(t)$
2. Integrate $\tilde\omega'(t)$ to get $\tilde\omega(t)$
3. Adjust $\tilde\omega'(t)$ and $\tilde\omega(t)$ to meet $\Omega$ constraint
4. Calculate $W(\tilde\omega(t)-\tilde\omega(0))$ and the new $\tilde\omega'(t)$
5. Repeat steps 2-4 until the solution converges

I am paraphrasing some of this and simplifying the notation a bit. Either way, I interpreted the above to mean start with some polynomial coefficients, fit the polynomial coefficients in $\tilde\omega'(t)$ to $\tilde\omega'(0)/W(\tilde\omega(t)-\tilde\omega(0))$ based on the coefficients you've got, adjust for constraints, and repeat. This does not appear to work and I don't really see how it would.

Am I interpreting this wrong? Is the algorithm in this paper referring to some well-known numerical differential equation solution approximation method?

Edit

$W(x)$ is a weighting function, specifically a sidelobe taper window. Assume real and symmetric. For example, a raised cosine window, $$W(x) = 1+\frac{1-\alpha}{\alpha}\cos(Cx)$$ with $C$ being some constant to make the span work out right. With $\alpha=25/46$ this is the ever-popular Hamming window.

Answer 1

The proposed scheme is a variation of a Picard iteration as described here.

However the proposed scheme deviates a little bit from a "standard" procedure especially when

1. You need to rescale the function
2. You need to approximate the solution to get back into the polinomial subspace.

Therefore to prove convergence we can directly apply the Banach fixed point theorem for a suitable fixed point operator and prove that this operator is a contraction.

Following the steps you proposed in the question and the comments we define

$$ \varphi(\nu)=\frac{\mu(\nu)\nu(0)}{W\left(\mu(\nu)\int_0^t\nu(s)ds\right)} \quad\text{with}\\ \mu(\nu)=\frac{\Omega}{\int_{-T/2}^{T/2}\nu(s)ds} $$

Here $\nu$ plays the role of $\omega'$ and the operator $\varphi$ covers steps 1.-4. without the approximation. ($\mu$ incorporates the rescaling into the operator).

I played around a bit and although lengthy it seems doable to prove that

$$ \max|\varphi(\nu_0)-\varphi(\nu_1)|\leq C(\max(W),\min(W),\lambda_W,T,\Omega)\cdot \max|\nu_0-\nu_1| $$

where $\lambda_W$ is the Lipschitz constant of $W$. (If you are interested in exact value of $C$ I found, I can post it here as well). This shows that for some choices of $\Omega,T,W$ the constant $C$ should be smaller than $1$ and the scheme does indeed converge.

Depending on your choice of polinomial approximation you will generally speaking get an additional constant $C_P$ that has to be taken into account namely

$$ \max|P(\varphi(\nu_0))-P(\varphi(\nu_1))|\leq C_P\cdot \max|\varphi(\nu_0)-\varphi(\nu_1)| \\ \leq C_P \cdot C(\max(W),\min(W),\lambda_W,T,\Omega)\cdot \max|\nu_0-\nu_1| $$