Finding the unknown in an ODE - numerical method

Question

Say I have the following ODE: $\frac{d^2x}{dt^2} = k \cdot f(x,t)$ where $k$ is an unknown constant and I'm trying to solve for $k$ knowing that at $x = x_f \rightarrow \frac{dx}{dt} = v_f$. Also the initial conditions are: $t = 0 \rightarrow x = 0 \,\, \& \,\, \frac{dx}{dt} = 0$.

Furthermore, I'm attempting to solve this numerically rather than analytically.

So far I've come up with the following:

Let $\frac{d^2x}{dt^2} = \frac{dv}{dt}$ therefore we can now apply Euler's method (bear in mind that $\frac{dx}{dt} = v$):

$v_{n+1} = v_n + k \cdot f(x,y) \cdot \Delta t$

$x_{n+1} = x_n + v_n \cdot \Delta t$

$t_{n+1} = t_n + \Delta t$

where $\Delta t \approx 0$

The initial conditions have already been stated in the first paragraph. Having defined the recurrence relationship we can now iterate and plot a graph of $x(t)$. If I look at the value of $v$ at $x = x_f$ I clearly get a different result depending on the value of $k$. So the point is to find that value for $k$ such that when $x = x_f \rightarrow v = v_f$.

The "simple" solution would be trial and error, i.e, to try different values for $k$. If you exceed $v_f$ decrease $k$, otherwise, increase it and keep doing so until you're satisfied with a good enough number.

This method is the one I'm trying to avoid and I'm seeking for something a lot quicker and smarter. Any ideas?

Thanks

Answer 1

You can adopt Newton's method to find the value of $k$. Basically you have a function $v_{x_f}(k)$ which is the final velocity at $x_f$. You want to solve the equation $v_{x_f}f(k)-v_f=0$. Recall that Newton's method states that $$k_{n+1}=k_n-\frac{v_{x_f}f(k_n)}{v'_{x_f}f(k_n)}$$ Depending on $f$, you might be able to find the derivative in the denominator exactly; otherwise, you will have to approximate the derivative given $k_n$ and $k_n+\delta$ for a small $\delta$.

Answer 2

You can also use the existing boundary value problem solvers. For the one in python scipy (Matlab has similar facilities) you would set

def ode(s, x, p): k,T = p; return [ x[1], T*k*f(s*T,x[0]) ]
def bc(x0, x1, p): k,T = p; return [ x0[0], x0[1], x1[0]-xf, x1[1]-T*vf ]

s = np.linspace(0,1, 11);
x = [ s*xf, s*vf ];

res = solve_bvp(ode, bc, s, x, p=[1.,1.])
print res.message

s = np.linspace(0,1, 302);
x = res.sol(s);
k, T = p;
print "k=",k," T=",T
plt.plot(s*T, x[0])

Depending on $f$, this may fail, then one has to try with a better initial guess for the curve $x$ and the parameters.

The variable end point of the original problem was transformed into a problem with fixed integration interval by introducing the parameter $T$ and considering the ODE for $u(s)=x(sT)$, that is, $u'(s)=Tx'(sT)=Tk\,f(sT,u(s))$.