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 /.
MapThread[Rule, {variables, ylist}];(*Equivalent to step*
f/.Thread[Rule[variables,ylist]]*)
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.
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? $\endgroup$ – Alex Trounev Sep 18 at 5:26Sin[(Pi z)/L]
orSin[(2 Pi z)/L]
in the second integral? $\endgroup$ – Alex Trounev Sep 18 at 8:05Sin[(Pi t)/L]
andSin[(2 Pi t)/L]
. $\endgroup$ – Alex Trounev Sep 19 at 4:09