# ODE with singularity

Here is a simple example of singular ODE.

NDSolve[{R'[s] == 1/s (-R[s] + T[s]),
T'[s] == 1/s (R[s] - T[s]) - 3 s R[s]^(-2),
R[0] == T[0], R[1] == 1}, {R, T}, {s, 0, 1}]


and not surprisingly this outputs singularity error messages. What I am confused about is that I tried changing the range from 0.0001 to 1 to avoid singularity at s=0, and it still outputs singularity errors. The solutions for T and R are known to be non-zero in this range. I wonder why this error happens, and I also wonder if there is a rule of thumb to numerically evaluate ODEs with this kind of singularity which I face oftentimes.

To solve the system of differential equations using NDSolve in the Wolfram Language, you can follow these steps:

sol = NDSolve[{
R'[s] == 1/s (-R[s] + T[s]),
T'[s] == 1/s (R[s] - T[s]) - 3 s R[s]^(-2),
R[0] == T[0],
R[1] == 1
}, {R, T}, {s, 0, 1}]


Here's what each part of the code does:

1. Differential Equations:

• The system of differential equations is specified using NDSolve.
• R'[s] == 1/s (-R[s] + T[s]) represents the derivative of R with respect to s.
• T'[s] == 1/s (R[s] - T[s]) - 3 s R[s]^(-2) represents the derivative of T with respect to s.
2. Initial Conditions:

• Initial conditions are provided for both R and T at s = 0.
• R[0] == T[0] specifies that R and T are equal at s = 0.
3. Boundary Conditions:

• A boundary condition is provided for R at s = 1.
• R[1] == 1 specifies that R equals 1 at s = 1.
4. NDSolve Function:

• NDSolve is the function used to solve the differential equations numerically.
• It takes the system of differential equations, initial conditions, and the range of s ({s, 0, 1}) as arguments.
• The result is stored in the variable sol.

After executing this code, sol will contain the solution to the system of differential equations over the specified range of s. You can then use Plot or other visualization functions to analyze and visualize the solution if needed. For example:

Plot[Evaluate[{R[s], T[s]} /. sol], {s, 0, 1}, PlotLegends -> {"R[s]", "T[s]"}]


This will plot the solutions for R[s] and T[s] over the range s = 0 to s = 1.

• Your code is identical with the OPs and does not address the question asked. Commented Apr 2 at 18:30
• Perhaps eps = 0.0001; NDSolve[{R'[s] == 1/s (-R[s] + T[s]), T'[s] == 1/s (R[s] - T[s]) - 3 s R[s]^(-2), R[eps] == T[eps], R[1] == 1}, {R, T}, {s, eps, 1}, Method -> {"Shooting", "StartingInitialConditions" -> {R[eps] == 1, T[eps] == 1}}]? Commented Apr 2 at 18:34
• @MichaelE2 Thanks! So the idea is to specify the method as "shooting"? Can I assume this strategy works in general for this kind of problem? Commented Apr 2 at 19:02
• Say, for higher order equations? Commented Apr 2 at 19:03
• @physgj When you have boundary conditions at different values of s, then you have what is called a "boundary-value problem" (BVP). Two methods of solution are available in NDSolve, "Shooting" and "Chasing". You can search for them in the help center for details on their options. NDSolve was using "Shooting" with automatically chosen starting values, but they led to errors and the method failed (like using Newton's method with a bad starting position, if you know about that). Anyway "Shooting" works with any order/dimension of system, but finding starting values can be hard. Commented Apr 2 at 19:52