I am trying to solve the ODE :

$\frac{d^2T}{dx^2} = \omega_1 T+\omega_2 T^2$ + using different numerical methods. I have tried the following discretizations so far and none of them seem to be grid independent:

Finite Difference Method

$\frac{T_{i+1}^k + T_{i-1}^k - 2T_i^{k+1}}{\Delta x^2} = \omega_1(T_i^2)^k+\omega_2T_i^{k+1}$

Another Finite Difference Method with relaxed non-linearity

$\frac{T_{i+1}^k + T_{i-1}^k - 2T_i^{k+1}}{\Delta x^2} = \omega_1T_i^kT_i^{k+1} + \omega_2T_i^{k+1}$

A finite volume discretization using Picard's method of accounting non-linearities

A general finite volume scheme the discretization of which I shall provide soon, however this also doesn't turnout to be grid independent.

The BCs for this problem are as follows: $\frac{dT}{dx}=0$ at $x=0$ and $T=1$ at $x=1$

Any insights on solving this using a consistent scheme? Thanks in advance.

Edit 1: I took the analytical results of this ODE from this paper and wrote a code to obtain solutions from the three discretizations. I don't find the solution to be constant no matter how small a grid resolution I go to. Ideally speaking it must attain constanty after a certain cut off making it consistent with the parent ODE. My guess is the non-linearity causes all the trouble.

Edit 2: So regarding the indices and the notation I am using, $k$ denotes the iteration number (I am solving the resulting system of equations using Jacobi's method) and $i$ refers to the spatial index. Regarding the code, you can have a look at the following code that implements the three different schemes I have implemented so far using MATLAB:

 %============================POROUS FIN SOLUTION==========================%
%%=============================NON-LINEAR FDM============================%%
for k=1:100
    for i=2:mesh-1
hold on;
%==============================LINEARIZED FDM============================%% 
mesh = 10
for k=1:100
    for i=2:mesh-1
        t(i) = (t_old(i+1) + t_old(i-1) - (w2*dx^2)*t_old(i))/(2 + t_old(i)*w1*dx^2);
% source term variables - add routine to compute based on other params
w1=1.1422; w2=0.5710;
% mesh: there will be (n-1) cells 
n = 6;
% mesh resolution
dx = 1 / (n-1);
% theta
t = zeros(1,n+1); 
% impose right boundary BC
t(end) = 1;
% outer iterations for the main loop
iter = 100; 
%%==========================THE JACOBI ITERATOR==========================%%
for k = 1 : iter
    t_old = t;
    for i = 2 : n
        sc = Sc(w1,w2,t_old(i));
        sp = dSdt(w1, w2, t_old(i));
        if (i == 2)
            t(i) = update_cell_1(t_old(i+1), sc, sp, dx);
        elseif (i==n)
            t(i) = update_end(t_old(i+1), t_old(i-1), sc, sp, dx);
            t(i) = update(t_old(i+1), t_old(i-1), sc, sp, dx);
    t(1)=t(2); %impose the BC at the left end
%%============================PLOTTING AND VIS===========================%%
%construct x-vector
xs = dx/2;
xe = 1 - dx/2;
x = linspace(xs, xe, n-1);
x = [0,x,1];
legend('FDM-NL', 'FDM-L','FVM-L')
hold off;
%%=============================LOOP UPDATES==============================%%
%main update
function adv = update(te, tw, sc, sp, dx)
   adv = (te + tw - sc * dx^2)/(2 + sp * dx^2); 
%cell 1 update (to account for BC)
function cell_1 = update_cell_1(te, sc, sp, dx)
    cell_1 = (te - sc * dx^2)/(1 + sp * dx^2);
%cell n update (accounts for the shrunk ghost-cell
function adv_end = update_end(te, tw, sc, sp, dx)
    adv_end = (2 * te + tw - sc * dx^2)/(3 + sp * dx^2);
%%==========================SOURCE TERM FUNCTIONS========================%%
%original source term
function s = S(w1, w2, t)
    s = w1 * t^2 + w2 * t; 
% derivative and the Sp term 
function sp = dSdt(w1, w2, t)
    sp = 2 * w1 * t + w2;
%constant term 
function sc = Sc(w1, w2, t)
    sc = S(w1, w2, t) - (dSdt(w1, w2, t)) * t;

Likely that I don't have converged results as I haven't set a convergence condition for the loop. However, also consider the results from the code represented below: Comparison of solutions This hints that the solution I might have obtained is correct, nonetheless a thorough convergence analysis needs to be done.

Edit 3: (Extremely sorry for all the missing pieces of data and a very incomplete question, had some time crunch!) Now to get to what exactly I am puzzled about is the change in the solution with the change in the grid spacing (for a given method). Let's say I use my finite difference code solver and change the grid resolution, here are some results: Compare plots]2 The FDM solution changes with small changes in grid resolution. Now that's fine with small incremental increase in the grid points. What I can't digest is what happens when I increase the mesh points to something large. For instance consider the solution when I put mesh as 100: 100 points Note: This is the case when I use the very first method mentioned in my question (for clarity, $i$ represents the spatial grid point, and $k$ represents the Jacobi iteration number, given we use the Jacobi method for a solution to the resulting linear system of equations)

Closing this with this final edit I observed in my code that other methods simply need a greater number of iterations to converge to the same solution. I pre-set the number of iterations to a fixed value instead of a while-like loop which is causing all this discrepancy. Thank you everyone for being patient and offering your pieces of advice and solution.

  • $\begingroup$ First, I think there is an issue with the subscripts and superscripts in your finite difference equations, at least for the left-hand sides. What do you mean by "grid independent" ? Your result will always depend on the grid refinement you use, up to a certain level of refinement beyond which the solution can be considered to have converged. This level depends on the schemes used and on the problem. $\endgroup$
    – Laurent90
    Mar 30, 2021 at 18:17
  • $\begingroup$ What are the indices $i$ and $k$? Is $i$ a spatial grid point index and $k$ the iteration of a nonlinear scheme? $\endgroup$ Mar 30, 2021 at 22:29
  • $\begingroup$ @WolfgangBangerth yes precisely. $\endgroup$ Mar 31, 2021 at 12:04
  • $\begingroup$ @Laurent90 you are right in thinking so. However based on this discretization of the problem in the manner I have presented here, the analytical solution (obtained using Adomain decomposition) is obtained only for a particular grid resolution. If I switch to any other grid resolution the solution almost collapses. Ideally finer meshes must result in constancy of the solution without much change, but in this case it changes everytime with no real grid independence. $\endgroup$ Mar 31, 2021 at 12:07
  • 1
    $\begingroup$ As for the lack of grid independence: Show us some of the solutions you get. Do your fixed point iterations converge on every mesh? How do you assess that they converge? If they do, how do solutions look like on different meshes? My best bet is that you have a bug in your code, but it is hard to tell without actually seeing anything. $\endgroup$ Mar 31, 2021 at 13:51

1 Answer 1


To be honest, I do not understand what you are asking about. I suspect, what you are concerned about is the fact that when you increase the density of your grid points, the solutions drastically change, insteaded of gradually improving and converging to the actual solution. So I could only offer my attempt to find a method for numerical solution that does not have such a deffect. I simply chose to implement your first method, the standard finite difference method, but not just because of its simplicity, but because it coincides with a method based on conceptual discretization of the equation, rather than a purely numerical approach.

First, the equation is the Euler-Lagrange equation of a Lagrangian function and as such the solution is the optimal function that satisfies the principle of least action for a Lagrangian action.

Define the Lagrangian function $$L\Big(\,T, \, \frac{dT}{dx} \,\Big) \, = \, \frac{1}{2}\left(\frac{dT}{dx}\right)^2 \, + \,\frac{\omega_2}{3}\,T^3 \, + \, \frac{\omega_1}{2}\, T^2 $$ Then, one can define the action $$S[T] \, = \, \int_{0}^1 \, L\left(\,T, \,\, \frac{dT}{dx} \,\right)dx$$ Then the principle of least action says that $$\delta S[T] = 0$$ when the Euler-Lagrange equations $$\frac{d}{dx}\left( \frac{\partial L}{\partial \dot{T}}\left(\,T, \,\, \frac{dT}{dx} \,\right)\right) \, = \, \frac{\partial L}{\partial {T}}\left(\,T, \,\, \frac{dT}{dx} \,\right)$$ which is exactly the equation $$\frac{d^2 T}{dx^2} \, = \, {\omega_2}\,T^2 \, + \, {\omega_1}\, T $$ Now, instead of discretizing the equation itself, I am going to discretize the action functional: I am going to assume that the interval is divided into $N-1$ equal intervals it consists of $N$ equally spaced points $\{ ih \, : \, i = 0,1,...,N-1 \} $ where $h = \frac{1}{N-1}$. Then, we are looking for a discrete function $$T = \big(\, T_0, \, T_1, ..., T_{N-1} \,\big)$$ which is technically an $N$ dimensional vector, that is a critical point of the action function $$S(T_0, ..., T_{N-1}) = \sum_{i=0}^{N-1} \, L\left(\,T_i, \, \frac{T_{i+1} - T_i}{h} \,\right) $$ which is really a multivariable calculus problem, formulated as a system of equations: $$\frac{\partial S}{\partial T_i}(T_0, ..., T_{N-1}) \, = \, 0 \,\, \text{ for } \,\, i = 1, ..., N-2$$ $$T_0 = T_1$$ $$T_{N-1} = 1$$

If you calculate the partial derivatives above you get only two consecutive terms of the sum, i.e. $$\text{for } \,\, i = 1, ..., N-2$$ $$\frac{\partial S}{\partial T_i}\, = \,\frac{\partial }{\partial T_i}L\left(\,T_{i-1}, \, \frac{T_{i} - T_{i-1}}{h} \,\right) \, + \,\frac{\partial }{\partial T_i}L\left(\,T_i, \, \frac{T_{i+1} - T_i}{h} \,\right) \,= \, 0 \,\, $$ $$T_0 = T_1$$ $$T_{N-1} = 1$$ Explicitly, the equations become $$\text{for } \,\, i = 1, ..., N-2$$ $$\frac{T_i - T_{i-1}}{h^2} \, + \, \frac{T_i - T_{i+1}}{h^2} \, + \, \omega_2\,T_i^2 \, +\, \omega_1\,T_i \, =\, 0$$ $$T_0 = T_1$$ $$T_{N-1} = 1$$ and after solving for $T_{i+1}$ $$\text{for } \,\, i = 1, ..., N-2$$ $$T_{i+1} \, = \, \big(\omega_2h^2)\,T_i^2 \, +\, \big(\omega_1h^2 + 2\big)\,T_i \, - \, T_{i-1} $$ $$T_0 = T_1$$ $$T_{N-1} = 1$$ Observe that if we knew $T_{i-1}$ and $T_{i}$ we could calculate $T_{i+1}$. The problem is that we do not know $T_0$ or $T_1$ but $T_{N-1} = 1$ instead. Because I was lazy to go on with smart ways about figuring out how to do this, I simply used the midpoint binary division of initial conditions in order to get the final value of $T$ to be $1$ (actually, to be very close to $1$ within a chosen accuracy). Here is some python code I wrote to check how this works and plotted a solution:

import matplotlib.pyplot as plt
import numpy as np

def integrate(T_initial, N_points, omega1, omega2):
    T = np.empty(N_points, dtype=float)
    T[0] = T_initial
    T[1] = T_initial
    h = 1/(N_points - 1)
    for i in range(1, N_points-1):
        T[i+1] = (omega2*h**2)*T[i]**2 + (omega1*h**2 + 2)*T[i] - T[i-1]
    return T

def solve_bvp(N, omega1, omega2, epsilon):
    T_down = 0
    T_up = 1
    while True:
        T_mid = (T_down + T_up)/2
        T = integrate(T_mid, N, omega1, omega2)
        T_final = T[-1]
        if abs(T_final - 1) < epsilon:
            return T
        elif T_final < 1:
            T_down = T_mid
        elif T_final > 1:
            T_up = T_mid
##### Parameters:  
omega1 = 1.1422
omega2 = 0.5710
N = 21
accuracy = 0.00001

##### boundary value problem solver:

T = solve_bvp(N, omega1, omega2, accuracy)

##### plots:

x = np.arange(0,N)/(N-1)
plt.plot(x, T)
for k in range(N):
    plt.plot(x[k], T[k], 'ro')

This is a plot with a grid of $21$ points:

enter image description here

and this is a plot with a grid of $51$ points:

enter image description here

So the solutions are very much alike with no major differences.

  • $\begingroup$ I tested this solution out for different mesh spacings and it works! Thanks for this solution. However, upon looking at your final discretization and mine, I don't spot any difference. The only point of difference between your code and mine is the way the resulting system of equations are solved. While I use the $T_i$'s as my subject you use the $T_{i+1}$'s and iterate until the end point reached the specified boundary value. $\endgroup$ Apr 2, 2021 at 16:34
  • $\begingroup$ @RSuryaNarayan Yes, indeed, as I pointed out in my post, I arrived at and implemented the same finite difference method as yours. However, I derived it not by simlpy discretizing the second derivative, but by discretizing an action functional, based on the Lagrangian structure of the equation. This actually has some extra good properties, like preservation of the symplectic structure of the phase space (position-momentum space). The other thing was to use a version of the Poincare map in order to adjust the solution to the boundary conditions. $\endgroup$ Apr 2, 2021 at 17:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.