Numerical methods (for ODE/PDE) that could take approximate solutions/good initial guesses, and further refine it to an certain accuracy

Question

I am currently playing with an old analog computer, which could solve time-dependent ODE/PDEs pretty fast, without time-stepping; thus there is no convergence issues caused by time-stepping because of its computing nature. But the problem with analog computer's solutions is that they are not accurate due to physical limitations. I am very curious that: is there any numerical methods/solvers which can take analog computer's approximate solution (over the time domain) to further process it, and generate a more accurate solution??

Let me give an example of solving second order ODE describing the motion a mass-spring damper. The equation is the following: $$ x'' = -0.2\cdot x' - 0.4\cdot x;\quad x(0)=1, x'(0) =0;\quad t_{stop} = 60s. $$ To solve the above equation on an analog computer, we need to map the above equation to an electrical system. Usually an analog computer could perform several basic arithmetic operation in the continuous-time domain, e.g. addition, subtraction, multiplication, integration etc. The output of an integrator represent an state-variable of the ODE; the input of that integrator represent the corresponding first-order time derivative. By configuring the basic computing blocks in feedback loops, we could map the equation as the following: (I use Simulink)

After you load the initial conditions onto the integrators, you can let the analog computer run and solve. If you measure the electrical signal at the output of integrator1, you will get the solution of $x(t)$ over the time domain:

But, due to the physical limitations (e.g. electrical noise, offsets), the solution of $x(t)$ is not accurate. What I am looking for is a numerical method that can take the above solution of $x(t)$ by analog computer, e.g. the solutions $x(t=1s), x(t=2s), x(t=3s), x(t=4s)... x(t=60s)$, start from these approximate solution points and further refine these solution $x(t=1s), ... x(t=60s)$ to a much higher accuracy.

(This second order ODE is just a simple case for illustration purpose; it happens to have analytic expression of solutions. The more general case would be nonlinear ODEs with no analytic solution.)

Thanks in advance!! Any thoughts and suggestions are greatly welcome and appreciated!!

Answer 1

If you have a good initial estimate, Newton's method is hard to beat. Quadratic convergence means that the number of accurate decimal (binary) places doubles with each iteration. This assumes that the first derivative is changing slowly between your estimate and the real solution, which means the second derivative times your error (between the estimate and the real answer) is small compared to the first derivative. From physical arguments you know your solution is a damped sine wave, so fit it to $A \cos (\omega t) \exp(-\lambda t)$ What you really need for Newton's method is estimates of $A, \omega, \lambda$, not estimates of $y(t)$ which is what you get from your circuit. $A$ is easy, it is $y(0)$. I would take $\omega$ from the last zero crossing I could identify easily and $\lambda$ from the ratio of the first peak to the starting amplitude.

Answer 2

In almost the same spirit as Ross Millikan's answer, let me assume that you know the initial conditions $x(0)=1$, $x'(0)=0$ and that the model is something like $$x(t)=e^{-\alpha t} \big(A \cos (\beta t)+B \sin (\beta t)\big)$$ and from data points $(t_i,x_i)$ taken from the analog computer, you want to refine the values of the four parameters $A,B,\alpha,\beta$ which appear in the model.

This can be considered as a nonlinear regression problem and the key issue is to get "reasonable" estimates.

The first condition gives $A=1$ and this is a definitive value (it is no more a parameter to be tuned).

Now, looking at the derivative $$x'(t)=e^{-\alpha t} (\cos (\beta t) (\beta B-\alpha A)-\sin (\beta t) (A \beta +\alpha B))$$ gives $\beta B-\alpha A=0$ which means that we can impose $B=\frac{\alpha} {\beta}$.

All of the above makes the model to be $$x(t)=e^{-\alpha t} \left(\frac{\alpha \sin (\beta t)}{\beta }+\cos (\beta t)\right)$$ and we are just left with two parameters $\alpha,\beta$ for which we need estimates. Looking at $$x'(t)=-\frac{\left(\alpha ^2+\beta ^2\right) e^{-\alpha t} \sin (\beta t)}{\beta }$$ we see that the first minimum of $x(t)$ will correspond to $\beta t_*=\pi$; this gives an estimate for $\beta$. At this point, we have $$x(t_*)=-e^{-\frac{\pi \alpha }{\beta }}=-e^{-\alpha t_*}$$ from which the estimate of $\alpha$ is easily extracted.

Now, we have all the required elements to start the nonlinear least square fit of the data.

Just looking at the plot in the post, using the fact that the first minimum corresponds more or less to $t_*=5$, $x=-0.6$, we obtain as estimates $\beta_0 = \frac{\pi }{5} \approx 0.628319$ and $\alpha_0= \frac{1}{5} \log \left(\frac{5}{3}\right)\approx 0.102165$ while the exact values should be $\beta=\frac{\sqrt{39}}{10}\approx 0.624500$ and $\alpha=\frac 1 {10}$. The nonlinear regression converges in a couple of iterations.