2
$\begingroup$

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.

enter image description here

(*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 /. 
      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);
    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]
|improve this question
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Alex Trounev Sep 18 at 5:26
  • $\begingroup$ 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. $\endgroup$ – mozeq Sep 18 at 7:58
  • $\begingroup$ Should there be Sin[(Pi z)/L] or Sin[(2 Pi z)/L] in the second integral? $\endgroup$ – Alex Trounev Sep 18 at 8:05
  • $\begingroup$ No, they're not in the integral. I've updated the figure of the equation. $\endgroup$ – mozeq Sep 18 at 8:39
  • $\begingroup$ 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]. $\endgroup$ – Alex Trounev Sep 19 at 4:09
2
$\begingroup$

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 enter image description here this processing allows it to be integrated.

enter image description here 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 enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Alex Trounev Sep 19 at 10:30
  • $\begingroup$ This makes the answer even more perfect $\endgroup$ – mozeq Sep 19 at 12:23
  • $\begingroup$ We can still improve the code, make it faster. I compared the three codes for speed - see update to my answer. $\endgroup$ – Alex Trounev Sep 19 at 14:45
2
$\begingroup$

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]

Figure 1 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. Figure 2

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}]]

Figure 3

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.

|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – mozeq Sep 19 at 7:17
  • $\begingroup$ OK! Do you have an exact solution for comparison? $\endgroup$ – Alex Trounev Sep 19 at 7:32
  • $\begingroup$ Actually it's not necessary to add ?NumberQ, without it NDSolve will be much faster. $\endgroup$ – xzczd Sep 19 at 8:53
  • $\begingroup$ 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 $\endgroup$ – mozeq Sep 19 at 9:57
  • 1
    $\begingroup$ @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. $\endgroup$ – xzczd Sep 21 at 7:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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