How to solve a non linear ODE with Newton's method?

Question

Im trying to solve this ODE using the Newton method (for non-linear equations using Jacobian Matrix): $$ \frac{u''(x)}{(1+u'(x)^2)^{3/2}}=\frac TK u(x)+\frac{wx(x-L)}{2K} $$ original ode image

$u(0) = u(L) = 0$ (let's take $L=9)$

$T=500$

$K= 5 \times 10^9$ $w = 100$

$u''$ and $u'(x)$ are defined as follows:

I don't know what is f1 and what is f1. How could this be implemented in Python? Anyone can help here please?

Answer 1

You get an automatic Newton solver by using the BVP solver scipy.integrate.solve_bvp. You program your equation as usual

def eqn(x,u): return [ u[1], (1+u[1]**2)**1.5*(T*u[0]+0.5*w*x*(x-L))/K ]
def bc(u0,uL): return [ u0[0], uL[0] ]

x = np.linspace(0,L,21);
a = 0.1
u = [ a*x*(L-x), a*(L-2*x) ]
res = solve_bvp(eqn, bc, x, u, tol = 1e-12); # the solution is small, scale 1e-6, thus reduce from the default tol=1e-3
print res.message
if res.success:
    x = np.linspace(0,L,700)
    u = res.sol(x)
    plt.plot(x,u[0])

resulting in the plot

Answer 2

Given a second order boundary value problem $y''(x)=f(x,y(x),y'(x))$, $y(a)=y_a$, $y(b)=y_b$, you can linearize it at a given approximate solution $y$ as \begin{align} y''(x)+Δy''(x) &= f(x,y(x)+Δy(x),y'(x)+Δy'(x)) \\ &\approx f(x,y(x),y'(x))+\partial_yf(x,y(x),y'(x))Δy(x)+\partial_{y'}f(x,y(x),y'(x))Δy'(x) \end{align} Now insert the difference quotients for $y$ to produce the coefficients and use the difference quotients for $Δy$ to establish the usual linear system for $$ Δy''(x) + a(x)Δy'(x)+b(x)Δy(x)=c(x) $$ leading to $$ (1+\tfrac12ha_k)Δy_{k+1}-(2-h^2b_k)Δy_k+(1-\tfrac12ha_k)Δy_{k-1}=c_k \\~\\ a_k=-\partial_{y'}f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h}) \\ b_k=-\partial_{y}f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h}) \\ c_k=f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h}). $$

---

With the given problem we have $$ f(x,y,y')=\frac{Ty-\tfrac12wx(L-x)}K(1+y'^2)^{3/2} $$ so that $$ \partial_{y}f=\frac{T}K(1+y'^2)^{3/2},\\ \partial_{y'}f=\frac{Ty-\tfrac12wx(L-x)}K\cdot 3(1+y'^2)^{1/2}y' $$ which now needs to be assembled into a code establishing the linear system for the $Δy_k$, which then is iteratively solved and used to update $y_k$.