# NDSolve for a system of nonlinear coupled pde

I try to solve a system($$n$$) of coupled nonlinear 1-d PDE's which is basically coupled sine-Gordon equations (through 2nd spatial derivative) in time and space. I am interested in the dynamics when I increase the force gradually. So I use the final answer of smaller force $$\beta$$ and use it as initial condition for a larger force $$\beta$$. I wrote a code for smaller force and it produces the result in a decent amount of time (depending on the $$n$$). I use the final result as initial condition in the second part of my code. Since both parts of my code contain basically the same type of equations and only different in initial conditions and force magnitude, I expect it to run the same way as the first part. But, unfortunately no matter what I tried, the second part takes very long time and runs forever. Here is my code:

n = 3; η = 0.1; β = 0.4; L = 10; T = 100; τ = T/10;
af = Array[f, n];
txaf[x_, t_] := Through[af[x, t]];
eqs = Join[
Array[D[f[#][x, t], t, t] + η D[f[#][x, t], t] -
D[Sum[i f[i][x, t], {i, 1, #}] + #*Sum[f[j][x, t], {j, # + 1, n}], x, x] +
Sin[f[#][x, t]] - β Tanh[t/τ] == 0 &,n]];
init = {Union[Thread[Drop[txaf[x, 0], {(n + 1)/2}] ==
0], {txaf[x, 0][[(n + 1)/2]] == 4 ArcTan[Exp[x]]}],
bc = {Thread[Derivative[1, 0][txaf][-L, t] == 0],
sol = NDSolve[{eqs, init, bc}, af, {x, -L, L}, {t, 0, T},
MaxStepSize -> 0.05]


using as intial condition:

f0[x_] := Thread[First[Evaluate[txaf[x, T] /. sol]]];
df0[x_] := Thread[First[Evaluate[Derivative[0, 1][txaf][x, T] /. sol]]];


second part of my code:

Clear[af, eqs, init, bc]
β2 = 0.5;
af = Array[f, n];
txaf[x_, t_] =.;
txaf[x_, t_] := Through[af[x, t]];
eqs = Join[
Array[D[f[#][x, t], t, t] + η D[f[#][x, t], t] -
D[Sum[i f[i][x, t], {i, 1, #}] + #*Sum[f[j][x, t], {j, # + 1, n}], x, x] +
Sin[f[#][x, t]] - β - (β2 - β)*
Tanh[t/τ] == 0 &, n]];
init = {Thread[txaf[x, 0] == f0[x]],
bc = {Thread[Derivative[1, 0][txaf][-L, t] == 0],
sol2 = NDSolve[{eqs, init, bc}, af, {x, -L, L}, {t, 0, T},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 1000, "MinPoints" -> 500, "DifferenceOrder" -> "Pseudospectral"}}]


you can copy and paste the code and just hit the run. In the mean time if you have any comment to make the code more efficient for higher number of $$n$$, it would be greatly appreciated.

• @bbgodfrey equation-solving tag should not be used in questions which don't deal with Reduce, Solve, NSolve, FindRoot, RSolve, LinearSolve et consortes, regardles of the fact that e.g. DSolve may use Solve behind the scene. Such a way of tagging is counterproductive. it produces more confusion than any potential advantages. Differential equations are too rich, they should not be confused with equation-solving. Nov 8, 2020 at 0:53
• @Artes The tag wiki explicitly includes NDSolve. Nov 8, 2020 at 1:26
• @bbgodfrey That tag wiki is simply harmful. Nov 8, 2020 at 2:00

As noted in the question, the first code segment does run to completion. However, NDSolve also issues two warning messages.

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

NDSolve::eerr: Warning: scaled local spatial error estimate of 37.88105237757757 at t = 100. in the direction of independent variable x is much greater than the prescribed error tolerance... >>

The first of these is due to the fact that init/.x -> L (or -L) does not match bc/.t -> 0. Correcting this should be straightforward. How serious the second warning is can be estimated from a plot of the solution.

Plot3D[Evaluate[Through[af[x, t]] /. Flatten@sol], {x, -L, L}, {t, 0, T},
AxesLabel -> {x, t, f}]


which does not show obvious signs of numerical instability but, nonetheless, may not be as accurate as desired.

Not surprisingly, the second calculation also flags initial and boundary conditions as inconsistent. Additionally, it warns,

NDSolve::mxsst: Using maximum number of grid points 1000 allowed by the MaxPoints or MinStepSize options for independent variable x. >>

As noted in the question, the calculation is exceedingly slow. In fact, several minutes are required to reach t == 1. This is due to the Method choice, "DifferenceOrder" -> "Pseudospectral", which works best for periodic boundary conditions. Simply deleting this option greatly decreases runtime, but still requires a great deal of memory. To overcome this latter problem, one can reduce MaxPoints to 500, so that NDSolve becomes

sol2 = NDSolve[{eqs, init, bc}, af, {x, -L, L}, {t, 0, 100}, Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 500}}]


which, after several minutes produces an answer similar to that from the first computation.

Plot3D[Evaluate[Through[af[x, t]] /. Flatten@sol2], {x, -L, L}, {t, 0, T},
AxesLabel -> {x, t, f}]


In summary, modifications to the options used by NDSolve lead to results in a reasonable amount of time. However, as noted above, the initial and boundary conditions should be modified so that init and bc are consistent at {L, 0} and {-L, 0). This, in itself, may improve the accuracy computations sufficiently to eliminate the second warnings.

• Thanks a lot. It was very helpful. But, I still do not have a clear image that why the second part should be any different in terms of computational time than the first part nevertheless they use the same set of equations. Apr 6, 2016 at 14:45
• @AhmadSheikhzada The second code as written was much slower, because "DifferenceOrder" -> "Pseudospectral" was inappropriate for this computation. Even with that option deleted, the initial conditions for the second code were less smooth than for the first, so NDSolve used finer resolution, which increased run time by a factor of between 2.5 and 2.5^2. Worse, the correspondingly greater memory usage required results to be swapped to disk, which caused the computation to come to almost a halt. I recommend you focus now on eliminating the first warning message. Apr 6, 2016 at 14:55
• Thanks again. I removed the first part of errors and even the second error in the second part. I am wondering if this code is parallizable, because the ultimate goal is to do this for larger number of equations. Apr 6, 2016 at 15:03
• @AhmadSheikhzada Probably, but a lot of work may be involved to accomplish it, for instance by converting the PDE into a set of coupled ODEs. It all depends on the details of the equations. Apr 6, 2016 at 15:45