How to numerically set up to solve this differential equation?

Question

I have a 1-d differential equation: $$\frac{\mathrm{d}f}{\mathrm{d}\theta} = c(\mathrm{max}(\sin\theta,0)-f^4)~.$$ I am given periodic boundary condition, i.e. $f(\theta) = f(2\pi+\theta)$. How would I set up a discretised form of this equation to solve for $f(\theta)$?

Answer 1

This problem is harder than a traditional finite-differencing problem. Traditional "time"-stepping methods will not work because your boundary conditions are not in the form of an initial-value problem.

The way you should think about this is setting up a nonlinear equation then performing Newton's method. From your periodicity requirement, you only need to consider points in $[0,2\pi)$. Generate an equispaced grid with $N+1$ points $\theta_i, \ i=0,\dots,N$ and $N$ unknowns $f_i, \ i=0,\dots,N$. Let $f$ be the vector of unknowns and let $h$ be the step size between the grid points. We can then write our problem as $$G(f) = 0,$$ where $$G(f)_i= \frac{f_{i+1}-f_{i-1}}{2h}-c(\max\{\sin\theta_i,0\}-f_i^4).$$ In your implementation, make sure that the equations for $i=0$ and $i=N$ are properly adjusted to account for the periodic conditions. Applying Newton's method with a forward difference Jacobian-vector product and a Krylov solver of your choice should do the trick for $N$ in the thousands or so.

This is all assuming that a solution exists and that your initial guess is close enough. It is also not obvious that the continuous problem has a solution, so it is possible that as you take $h\to0$, things stop making sense since you are trying to find an answer that isn't there.

Answer 2

Any higher-order numerical method will experience a singularity at $θ=0$ and $θ=\pi$, as the error estimator in the step size control depends on the smoothness of the derivatives of the right side up to some order connected to the order of the method. A kink like here will be "seen" as a large oscillation in these derivatives, necessitating emergency procedures.

Thus it is best to separate the two cases and integrate them separately, taking the initial values of the second piece from the corresponding values of the first segment.

---

In this case you have the switch between the two modes at equally spaced points, so you can fold back the second half interval and consider $$F(θ)=[f(θ),f(2\pi-θ)]$$ with $\dot F=[c(\sinθ-F_1^4), cF_2^4]$ for $θ\in [0,\pi]$ with the periodicity/boundary condition $F_2(0)-F_1(0)=0$, $F_2(\pi)-F_1(\pi)=0$.

Use that to set up your preferred boundary value solver

from scipy.integrate import solve_bvp;

c = 0.5;
def F_ode(t,F): return [ c*(np.sin(t)-F[0]**4), c*F[1]**4 ]
def F_bc(y0,y1): return [ y0[1]-y0[0], y1[1]-y1[0] ]

x = np.linspace(0,np.pi, 11)
F = [ 1+0*x, 1+0*x ]

res = solve_bvp(F_ode, F_bc, x, F, tol=1e-9)
print res.message
if res.success:
    F=res.sol(x);
    for xk, Fk in zip(x,F.T): print "%15.10f -> %15.10f  | %15.10f -> %15.10f"%(xk,Fk[0],2*np.pi-xk, Fk[1])  
    x = np.linspace(0,np.pi, 150);
    F = res.sol(x);
    plt.plot(x,F[0],2*np.pi, F[1])
    plt.grid(); plt.show()

which finishes successfully and produces the plot

and function table $x\to f(x)$:

   0.0000000000 ->    0.5366355475  |    6.2831853072 ->    0.5366355475
   0.3141592654 ->    0.5479126886  |    5.9690260418 ->    0.5503329371
   0.6283185307 ->    0.6020560795  |    5.6548667765 ->    0.5655453511
   0.9424777961 ->    0.6857306918  |    5.3407075111 ->    0.5825927129
   1.2566370614 ->    0.7794739956  |    5.0265482457 ->    0.6019011175
   1.5707963268 ->    0.8617720178  |    4.7123889804 ->    0.6240537356
   1.8849555922 ->    0.9166873730  |    4.3982297150 ->    0.6498755212
   2.1991148575 ->    0.9383800806  |    4.0840704497 ->    0.6805822042
   2.5132741229 ->    0.9284736894  |    3.7699111843 ->    0.7180612125
   2.8274333882 ->    0.8908228729  |    3.4557519189 ->    0.7654512168
   3.1415926536 ->    0.8284926309  |    3.1415926536 ->    0.8284926309

For the more detailed case description in question How to remove this numerical artifact? which gives $c=33.33$ here, you get the solution