in fluid mechanics I came across with the so-called Stream function that is capable of describing waves and its underlying kinematics quite accurately. However, its implementation based upon Newton-Raphson is tedious.
Short outline of the idea (see Dean 1965, Dalrymple 1974): The objective function, OF, is now to be minimised with respect to X and L (i.e. a first order Fourrier coefficient X and the wavelength L) . $$ OF = 1/L* \int_0^L (Qi-Q)^2dx \quad\quad(1) $$ with Qi the Bernoulli quation evaluated on each point on the wave profile $$ Qi=\eta + 1/(2*g)*((u-c)^2+w^2) - c^2/(2*g) $$ and Q its average $$ Q=1/L * \int_0^L Qi dx $$ where g=9.81, c=L/T the wave speed, and the velocity in x-direction (horizontal) $$ u=-X*exp(2*\pi/L*\eta)*cos(x) $$ and the velocity in z-direction (vertical) $$ w=-X*exp(2*\pi/L*\eta)*sin(x) $$ x is the respective value on the x-axis and ranges from 0 to 2*pi (one wavelength). The surface profile is given by $$ \eta = (1/c*\Psi - 1/(2*\pi/L) * (X/(2*\pi/L)*exp((2*\pi/L)*\eta).*cos(x))) $$ and is a vector of length x. The constant (and very small) Streamfunction parameter is $$ \Psi = 1/L* \int_0^L X/(2*\pi/L) * exp((2*\pi/L)*\eta) .* cos(x) dx $$
Now, the solution procedure is, given an initial guess for X, L and $\eta$, which is close to the desired solution, minimise equation (1) by changing X and L (not $\eta$). So, $$ dOFdX = 1/L * \int_0^L (Qi-Q)*(dQidX - dQdX) dx $$ $$ dOFdL = 1/L * \int_0^L (Qi-Q)*(dQidL - dQdL)dx $$ where the derivatives need to be determined using the chain rule: $$ dQidX = dQideta.*detadX + dQidu.*(dudX + dudeta.*detadX) + dQidw.*(dwdX + dwdeta.*detadX) $$ $$ dQidL=dQideta.*detadL + dQidu.*(dudL + dudeta.*detadL) + dQidw.*(dwdL + dwdeta.*detadL) + dQidc.*dcdL; $$ and dQdX, dQdL the respective averages over the wavelength L.
Now, the routine looks like this:
1) find an initial guess which is close to the solution. $$ \eta = -0.01*cos(x) $$ where 0.01 is the wave amplitude in metres. X and L are 10% off the exact solution.
2) determine the X and L for the next iteration using Newton-Raphson (and a damping coefficient $\alpha$), so $$ X'= -\frac{d^2OF}{dP^2} * \frac{dOF}{dX} $$ $$ X^{j+1} = X + \alpha * X' $$ (j: iteration step) $$ L' = -\frac{d^2OF}{dP^2} * \frac{dOF}{dL} $$ $$ L^{j+1} = L + \alpha * L' $$ where $$ \frac{d^2OF}{dP^2} $$ is the Hessian matrix, i.e. $$ \begin{pmatrix} \frac{d^2OF}{dX^2} & \frac{d^2OF}{dXdL} \\ \frac{d^2OF}{dLdX} & \frac{d^2OF}{dL^2} \end{pmatrix} $$
3) calculate the surface profile $\eta$ based on the new X and L value.
4) calculate the $\Psi$ value based on the new $\eta$ value and new X and L.
5) go to 2) until OF is sufficiently small.
Problem: $\eta$ is an unkown itself. Minimising OF will reduce the wave amplitude as this reduces OF.
Observations so far: The convergence of the wavelength L works very well. But the convergence of X which also determines the wave amplitude goes towards the wrong value, which is always smaller than the correct value (although the "convergence pattern" looks perfectly fine). Also the gradient after the first iteration has the wrong sign (should be positive, is negative). Changing the sign does not help though. Consequently, the wave amplitude is too small. So, in the next iteration X would be reduced even further. Hence, there seems to be a coupling problem between eta and X. Providing a "correction function" that adjusts the wave height after each iteration makes things even worse. Even if I provided the correct surface profile and did not change it during the iteration, both values, X and L, would not be correct, so apparently a degree of freedom in eta must be retained. The gradient functions are analytically correct though.
Are there any ideas how to resume the issue of convergence. I have already checked the article "Condition under which newton raphson converges".
Any help is very much appreciated!
Best wishes and thanks, Tom.