# Trouble second-order ordinary differential equations with integral terms

I am a beginner of mma. During the learning of differential equations, I encountered a very difficult problem, a second-order system of ordinary differential equations with integral terms. This equation is very complicated to solve by using the built-in function of mma, so I plan to use the self-programmed RK4 method to solve it.After reading the RK4 function written by others, the code I wrote is as follows, but no matter what, I cannot get the result, which makes me very upset. I really hope that someone can help me to solve this problem. This is my first question, please forgive me for my unclear wording.

(*RK4*)
rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[{table, xlist, ylist, step, k1, k2, k3, k4}, xlist = tinit;
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;
xlist = tinit;
table = {{xlist, ylist}};
Table[k1 =
step*f /.
k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];
k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];
k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
xlist += step;
AppendTo[table, {xlist, ylist}];
{xlist, ylist}, nsteps];
table];

Here is the code I try

(*parameter*)
L = 100;
Cf[z_] := (p1*Sin[(Pi z)/L] + p2*Sin[(2 Pi z)/L])*RealAbs[p1*Sin[(Pi z)/L] + p2*Sin[(2 Pi z)/L]];

funclist = {p1, p2,
Sin[(Pi t)/L] - q1 - p1 -  NIntegrate[Cf[z]*Sin[(Pi z)/L], {z, 0, L}],
Sin[(2 Pi t)/L] - q2 - p2 - NIntegrate[Cf[z]*Sin[(Pi z)/L], {z, 0, L}]};
initials = {0, 0, 1, 1};
variables = {q1, q2, p1, p2};
init = 0;
final = 10;
nstep = 1000;
approx = rk4[funclist, variables, initials, init, final, nstep]
• In the code you have NIntegrate[Cf[z]*Sin[(Pi z)/L], {z, 0, L}] and in figure $\int_0^L Cf(z)\sin (\pi x/L)dz$ . Where is it right? Sep 18, 2019 at 5:26
• I‘m sorry, the "-Sin[Pi t/L]" isn't in the intergrate term. And The place where "x" is written in the figure should be "z".The code is right. Sep 18, 2019 at 7:58
• Should there be Sin[(Pi z)/L] or Sin[(2 Pi z)/L] in the second integral? Sep 18, 2019 at 8:05
• No, they're not in the integral. I've updated the figure of the equation. Sep 18, 2019 at 8:39
• You use RK4 module for an autonomous system. But in your system there is a clear dependence on time through Sin[(Pi t)/L] and Sin[(2 Pi t)/L]. Sep 19, 2019 at 4:09

In fact, Cf(z) is a nonlinear force, and q1 and q2 represent the first two modes of the vibrating system.I might want to solve for modes 5 - 10,That means there are 5-10 coupled equations.And Cf(z) is going to be 5 × 5 terms - 10 × 10 terms. For j=2 this processing allows it to be integrated.

But the result is very unfriendly, and although you can solve it with NDSolve, you can't solve it once the mode increases.The answer is the same as that of Alex Trounev. Thanks for his enthusiasm and answers

• This is a very interesting solution to the problem. As you can see, all solutions to the problem coincide. I have added an algorithm using NIntegrate[] to my answer Sep 19, 2019 at 10:30
• This makes the answer even more perfect Sep 19, 2019 at 12:23
• We can still improve the code, make it faster. I compared the three codes for speed - see update to my answer. Sep 19, 2019 at 14:45

In this case, it makes no sense to use RK4, you can use the standard solver and Gauss quadrature formulas for calculating integrals.

Get["NumericalDifferentialEquationAnalysis`"]; L = 100;
np = 60; points = weights = Table[Null, {np}]; Do[
points[[i]] = GaussianQuadratureWeights[np, 0, L][[i, 1]], {i, 1, np}]
Do[weights[[i]] = GaussianQuadratureWeights[np, 0, L][[i, 2]], {i, 1,
np}]
GaussInt[f_, z_] :=
Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]

Cf[z_, p1_, p2_] := (p1*Sin[(Pi z)/L] + p2*Sin[(2 Pi z)/L])*
Abs[p1*Sin[(Pi z)/L] + p2*Sin[(2 Pi z)/L]];
F1[p1_?NumberQ, p2_?NumberQ] :=
GaussInt[Cf[z, p1, p2]*Sin[(Pi z)/L], z];
F2[p1_?NumberQ, p2_?NumberQ] :=
GaussInt[Cf[z, p1, p2]*Sin[(2 Pi z)/L], z]
eq = {q1'[t] == p1[t], q2'[t] == p2[t],
p1'[t] == Sin[(Pi t)/L] - q1[t] - p1[t] - F1[p1[t], p2[t]],
p2'[t] == Sin[(2 Pi t)/L] - q2[t] - p2[t] - F2[p1[t], p2[t]]};
ic = {q1[0] == 0, q2[0] == 0, p1[0] == 1, p2[0] == 1};
var = {q1, q2, p1, p2};

sol = NDSolveValue[{eq, ic}, var, {t, 0, 10}];
Plot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 10},
PlotLegends -> {"q1", "q2"}, AxesLabel -> Automatic]

We will explain how to use rk4 to solve this problem. Since the system of equations is not autonomous, we include t in the number of variables. We remodel rk4 a bit and compare the result with a standard solver.

(*RK4*) rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] :=
Module[{table, ylist, step, k1, k2, k3, k4},
step = N[(tfinal - tinit)/(nsteps)];
ylist = valtinit;

table = {ylist};
Table[k1 =
step*f /.
k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];
k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];
k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];
ylist += 1/6 (k1 + 2 (k2 + k3) + k4);

AppendTo[table, ylist];
ylist, nsteps];
table];

(*parameter*)L = 100;
funclist = {p1, p2, Sin[(Pi t)/L] - q1 - p1 - F1[p1, p2],
Sin[(2 Pi t)/L] - q2 - p2 - F2[p1, p2], 1};
initials = {0, 0, 1, 1, 0};
variables = {q1, q2, p1, p2, t};
init = 0;
final = 10;

nstep = 1000;
approx = rk4[funclist, variables, initials, init, final, nstep];
q1s = Table[{approx[[i, 5]], approx[[i, 1]]}, {i, 10, Length[approx],
20}];
q2s = Table[{approx[[i, 5]], approx[[i, 2]]}, {i, 10, Length[approx],
20}];
Show[Plot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 10},
PlotLegends -> {"q1", "q2"}, AxesLabel -> Automatic],
ListPlot[{q1s, q2s}]]

We see here a rather exact coincidence of numerical solutions.

We now show how to use NIntegrate[] to solve this problem.

funclist1 = {p1, p2, Sin[(Pi t)/L] - q1 - p1 - F11[p1, p2],
Sin[(2 Pi t)/L] - q2 - p2 - F21[p1, p2], 1};
F11[p1_?NumberQ, p2_?NumberQ] :=
NIntegrate[Cf[z, p1, p2]*Sin[(Pi z)/L], {z, 0, L},
Method -> {"Trapezoidal", "SymbolicProcessing" -> 0}];
F21[p1_?NumberQ, p2_?NumberQ] :=
NIntegrate[Cf[z, p1, p2]*Sin[(2 Pi z)/L], {z, 0, L},
Method -> {"Trapezoidal", "SymbolicProcessing" -> 0}];

approx = rk4[funclist1, variables, initials, init, final, nstep];

q11s = Table[{approx[[i, 5]], approx[[i, 1]]}, {i, 10, Length[approx],
20}];
q21s = Table[{approx[[i, 5]], approx[[i, 2]]}, {i, 10, Length[approx],
20}];

We see again here a rather exact coincidence of numerical solutions.

Show[Plot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 10},
PlotLegends -> {"q1", "q2"}, AxesLabel -> Automatic],
ListPlot[{q11s, q21s}]]

Compare the three codes for speed: 1) a code using rk4 and GaussianQuadratureWeights - 7.5 seconds; 2) code using GaussianQuadratureWeights and NDSolve - 1.038 sec; 3) code using rk4 and NIntegrate[] - 9.36 sec. If we remove ?NumberQ in Definition F1 and F2, as suggested by @xzczd , then the speed increases significantly: 1) - 3.2 s; 2) - 0.23 s.

• You hit the nail on the head.If you do the integration well, it's not hard to solve.But I can't think of a way to do that.Finally, thank you very much. Sep 19, 2019 at 7:17
• OK! Do you have an exact solution for comparison? Sep 19, 2019 at 7:32
• Actually it's not necessary to add ?NumberQ, without it NDSolve will be much faster. Sep 19, 2019 at 8:53
• Once we've done the integration, we don't have to do the reduction of the second order differential equation, we can just solve it.I think there is no exact solution to this equation Sep 19, 2019 at 9:57
• @CATrevillian Oops, I misunderstood your previous comment. The code becomes faster without ?NumberQ in this case because F1[…] and F2[…] are both calculated symbolically, thus the Sum[…] inside GaussInt[…] is calculated only once in F1 and F2 respectively, while with ?NumericQ the summation needs to be calculated numerically every time F1/F2 is called. Sep 21, 2019 at 7:34