Solve Second-Order ODEs Numerically

I have a complicated second-order ODEs. af(x)=b2f′(x)2+bpa+cf′(x)−cxf′(x)+σ22x2f″(x)af(x)=\frac{b}{2}f'(x)^2+\frac{bp}{a+c}f'(x)-cxf'(x)+\frac{\sigma^2}{2}x^2 f”(x).

This ODE has no analytical solutions, thus I want to solve it numerically. For example, I assume ALL constants to be 1 so the problem becomes f(x)=12f′(x)2+12f′(x)−xf′(x)+12x2f″(x)f(x)=\frac{1}{2}f'(x)^2+\frac{1}{2}f'(x)-xf'(x)+\frac{1}{2}x^2 f”(x). And combining with my initial conditions f(1)=3/8f(1)=3/8 and f'(1)=-1/2f'(1)=-1/2, f”(1)=0.

Yes I have three initial values, but may be not all should be used? Also in my problem setting, intuitively \lim_{x\rightarrow\infty}f(x)=0\lim_{x\rightarrow\infty}f(x)=0.

I am wondering if there is any numerical methods I can use to solve it? Or could anyone provide some references I can read?

=================

If you use Matlab, you may want to look at N. Trefethen’s Matlab Toolbox “Chebfun”. It uses spectral collocation combined with Newton root-finding techniques to solve fairly complicated problems. Here is the nonlinear ODEs learning-by-example demos webpage: chebfun.org/examples/ode-nonlin/index.html
– Dmoreno
2 days ago

=================

1

=================

The standard method (without looking at the differential equation here in detail and coming up with something fancy) is to rewrite a second order system in a coupled first order system and then solve this coupled first order system with standard methods (Newton, Runge-Kutta,…). In order to do that lets define
g(x)\equiv f'(x).g(x)\equiv f'(x).
With that your system becomes a coupled system of first order differential equations:

\begin{align}
f'(x)&=g(x),\\
\frac{\sigma^2}{2}x^2 g'(x)&=a f(x)- \frac b 2 g(x)^2-\frac{b p}{a+c}g(x)+c x g(x).
\end{align}\begin{align}
f'(x)&=g(x),\\
\frac{\sigma^2}{2}x^2 g'(x)&=a f(x)- \frac b 2 g(x)^2-\frac{b p}{a+c}g(x)+c x g(x).
\end{align}
For x\neq0x\neq0 one can rewrite this as:

\begin{align}
f'(x)&=g(x),\\
g'(x)&=\frac{2}{\sigma^2x^2}\left(a f(x)- \frac b 2 g(x)^2-\frac{b p}{a+c}g(x)+c x g(x)\right).
\end{align}\begin{align}
f'(x)&=g(x),\\
g'(x)&=\frac{2}{\sigma^2x^2}\left(a f(x)- \frac b 2 g(x)^2-\frac{b p}{a+c}g(x)+c x g(x)\right).
\end{align}

Integrating this system numerically gets involved if you have only one initial condition at x=1x=1. One has to use a shooting method to solve it as a boundary value problem. For that one should probably change the variables to compactify the integration domain: something like x\Rightarrow 1/ux\Rightarrow 1/u so that one has to integrate the system from u_0=1u_0=1 up to u_\infty=0u_\infty=0. On that domain one can start the numerical integration at u_0u_0 with the given initial value and one arbitrary one: this gives a solution at u_\infty=0u_\infty=0. Now one can use a minimizer to optimize the second initial condition until one gets the right boundary value (f(u_\infty=0)=0).

Thanks a lot. My problem has been refined to some extent.
– Daniel
2 days ago