0
$\begingroup$

Recently, I am doing a course project, and the professor forces us to use Mathematica; however, I am new to Mathematica, and here are my problems.

Setup

I am required to, eventually, numerically solve a group of differential equations which look like the form: $$\frac{dX(t)}{dt}=P(X(t))+F(t) $$ Here, $X(t)$ is the coloumn vector $(x_1(t),x_2(t),x_3(t))^T$, $P(X(t))$ is a polynomial of $X(t)$ which contains terms like $x_i(t)x_j(t)$.

$F(t)$ is the key term for this project. $F(t)$ is a column vector takes form $(F_1[x_1(t)]\cdot x_1^2(t),\,F_2[x_2(t)]\cdot x_2^2(t),\,F_3[x_3(t)]\cdot x_3^2(t))$, where $F_i[x_i(t)]$ is a upper bound function of a complicated function.

Here, I will just present $F_1[x_1(t)]$, since $F_2$ and $F_3$ are similar to $F_1$. I think I can deal with $F_2$ and $F_3$ once I know the way to deal with $F_1$.

$$F_1[x_1(t)]:=\int_0^{\frac{x_1(t)}{3}}\frac{(1+y)(2y-1)^2(2y-3)^2\sin^3(\pi y)[\Gamma(y-1)]^2\Gamma(-2y)}{(y-2)\pi^3}dy $$ Here $\Gamma$ is the usual gamma function. So far, we are given that $F_1(x)$ has poles at $x=\frac{15}{2}+3n$. Similarly to $F_2$ and $F_3$. Here, we only care about the value before the first pole.

Main Problem

Now, here is my question. In order to plug $F_1[x_1(t)]$ in to the differential equation, I tried to numerically compute $F_1(x)$ and define $F_1(x)$ as the following code:

I1[x_] := (1 + x)*(2*x - 1)^2*(2*x - 3)^2*Sin[(\[Pi]*x)]^3*
   Gamma[(x - 1)]^2 Gamma[(-2)*x]/((x - 2)*\[Pi]^3)

F1[x_] := NIntegrate[I1[t]*Boole[t < 2.5 - Exp[-NF2]], {t, 0, x/3}]

Since we only care $0$ to $15/2$, I set a Boole Function take the value from $0$ to $2.5-\text{Exp(-NF2)}$ $[\text{Note: }15/2\div 3]$, where NF2 is some large number which is given by professor.

After I define the $F_1(x)$, I plug the whole $F_1[x_1(t)]$ into the differential equation, then I get a lot of warning. Can someone help me to fix the issue? Thanks.

Attached is the graph of $F_1(x)$ and the error for directly plug-in procedure. $\alpha(t)$ is $x(t)$ here.Graph of F1Error1Error2Error3

Code

Here is the code I am working on.

NF3 = 40; NF2 = 24; Y2 = 
 1/2; Y3 = 0; DR22 = 2; DR32 = 1; DR23 = 1; DR33 = 3; \[Alpha]y = 0;
b1 = 41/3 + 8/3*Y2^2*NF2*DR22*DR32;
b2 = 19/3 - 4*NF2*DR32/3;
b3 = 14 - 4*NF3*DR23/3;
c1 = 199/9 + 8/3*Y2^4*NF2*DR22*DR32;
c2 = 35/3 + 49*NF2*DR32/3;
c3 = -52 + 76*NF3*DR23/3;
d1 = 88/3 + 32/3*Y2^2*NF2*DR22*DR32;
d2 = 24 + 16/3*NF2*DR32;
d3 = 9 + 3*NF3*DR23;
e1 = 9 + 6*Y2^2*NF2*DR22*DR32;
e2 = 3 + 4*Y2^2*NF2*DR32;
e3 = 11/3 + 4*NF3*Y3^2*DR23;

NC2 = 2;
NC3 = 3;
I1[x_] := (1 + x)*(2*x - 1)^2*(2*x - 3)^2*Sin[(\[Pi]*x)]^3*
   Gamma[(x - 1)]^2 Gamma[(-2)*x]/((x - 2)*\[Pi]^3);
F1[x_] := NIntegrate[I1[t]*Boole[t < 2.5 - Exp[-NF2]], {t, 0, x/3}];
I22[x_] := (NC2^2 - 1)/NC2 + (20 - 43 x + 32 x^2 - 14 x^3 + 4 x^4)*
    NC2/(2*(2 x - 1)*(2 x - 3)*(1 - x^2));
I23[x_] := (NC3^2 - 1)/NC3 + (20 - 43 x + 32 x^2 - 14 x^3 + 4 x^4)*
    NC3/(2*(2 x - 1)*(2 x - 3)*(1 - x^2));
H12[x_] := -11/2*NC2 + 
   NIntegrate[I1[t]*I22[t]*Boole[t < 1 - Exp[-NF3]], {t, 0, x/3}];
H13[x_] := -11/2*NC3 + 
   NIntegrate[I1[t]*I22[t]*Boole[t < 1 - Exp[-NF3]], {t, 0, x/3}];
Plot[F1[t], {t, 0, 8}]

A1 = 4*\[Alpha]1[t]*NF2*Y2^2*DR22*DR32;
A2 = 2*\[Alpha]2[t]*NF2*DR32;
A3 = 2*\[Alpha]3[t]*NF3*DR23;

odeg1 = 2*A1*\[Alpha]1[t]*F1[A1]/(3*NF2) + 
    D[\[Alpha]1[t], 
     t] == (b1 + c1*\[Alpha]1[t] + d1*\[Alpha]3[t] + 
      e1*\[Alpha]2[t] - 17*\[Alpha]y/3)*\[Alpha]1[t]^2;
odeg2 = 2*A2*\[Alpha]2[t]*H12[A2]/(3*NF2) + 
    D[\[Alpha]2[t], 
     t] == (-b2 + c2*\[Alpha]2[t] + d2*\[Alpha]3[t] + 
      e2*\[Alpha]1[t] - 3*\[Alpha]y)*\[Alpha]2[t]^2;
odeg3 = 2*A3*\[Alpha]3[t]*H13[A3]/(3*NF3) + 
    D[\[Alpha]3[t], 
     t] == (-b3 + c3*\[Alpha]3[t] + d3*\[Alpha]2[t] + 
      e3*\[Alpha]1[t] - 4*\[Alpha]y)*\[Alpha]3[t]^2;
Ncorrect = 
  NDSolve[{odeg1, odeg2, 
    odeg3, \[Alpha]1[3] == 0.00084, \[Alpha]2[3] == 
     0.00256, \[Alpha]3[3] == 
     0.0061}, {\[Alpha]1, \[Alpha]2, \[Alpha]3}, {t, 0, 20}];
Plot[Evaluate[{\[Alpha]1[t], \[Alpha]2[t], \[Alpha]3[t]} /. 
   Ncorrect], {t, 0, 20}, PlotRange -> 1, AxesLabel -> {t, Strength}, 
 ImageSize -> Large]
$\endgroup$
  • 1
    $\begingroup$ Please post a copy-pastable code, not only an image of a code, so that people can play with it on their computers. $\endgroup$ – corey979 May 27 '18 at 17:12
  • $\begingroup$ @corey979 The only thing one needs is the code in the main problem. Other codes are irrelevant. $\endgroup$ – Hamio Jiang May 27 '18 at 17:14
  • 1
    $\begingroup$ You seem to not understand how this site work. Please take the tour and then edit your question accordingly. $\endgroup$ – corey979 May 27 '18 at 17:19
  • $\begingroup$ @corey979 Apparently not. Thanks for the link. $\endgroup$ – Hamio Jiang May 27 '18 at 17:24
  • $\begingroup$ Can't analyse in detail, but try F1[x_?NumericQ] := NIntegrate[I1[t]*Boole[t < 25/10 - Exp[-NF2]], {t, 0, x/3}, MaxRecursion -> 50, WorkingPrecision -> 20]; and H12[x_?NumericQ] := -11/2*NC2 + NIntegrate[I1[t]*I22[t]*Boole[t < 1 - Exp[-NF3]], {t, 0, x/3}, MaxRecursion -> 50, WorkingPrecision -> 20]; ,same for H13 and Ncorrect = NDSolve[Rationalize[{odeg1, odeg2, odeg3, \[Alpha]1[3] == 0.00084, \[Alpha]2[3] == 0.00256, \[Alpha]3[3] == 0.0061}, 0], {\[Alpha]1, \[Alpha]2, \[Alpha]3}, {t, 0, 20}, MaxSteps -> 200, WorkingPrecision -> 20] $\endgroup$ – Akku14 May 27 '18 at 18:24
2
$\begingroup$

The use of x_?NumericQ which means the function evaluates only for numerical values should fix your problem. However, to compute F1, instead of computing the (indefinite) numerical integral one can turn it into a first-order DE:

$f(x)=\int_0^x g(t) \, dt \implies f'(x)=g(x)$ with $f(0)=0$

and solve this using NDSolve instead, which is much faster as it computes the values of F1 over the entire range in one go. Also, the solution is returned as an InterpolatingFunction so the issue that you faced disappears. But, because of the singularities of $g(x)$ at $x=0$ and $x=5/2$, one has to deal with these separately. To do this, perform series expansion about $x=5/2$ and subtract of the resulting pole there, (which leads to a logarithmic singularity).

|improve this answer|||||
$\endgroup$

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.