How to solve this differential equation with integral?

Question

How to solve this differential equation with integral?

$\omega = 1.2;$

$i1 = 25640;$

$ia = 4*10^6;$

$c2 = 1.4513*10^5;$

$t5 = 16000*\cos({\omega*t});$

$tpto = 2.8*10^6*{\theta'}(t);$

$cr(\omega)=\frac{1.8*10^7}{0.5}{\omega}^2\exp({-\omega})$

$R_{55}(t)=\frac{2}{\pi}{\int}_{-\infty}^{10}cr(\omega)\cos({\omega}t)d{\omega}$

$(i1+ia)*{\theta}^{''}(t)+{\int}_{-\infty}^{t}R_{55}(t-\tau){\theta}^{'}(\tau)d{\tau} +(c2)*{\theta}(t)=t5+tpto$

$(i1+ia)*{\theta}^{''}(t)+{\int}_{0}^{t}R_{55}(t-\tau){\theta}^{'}(\tau)d{\tau} +(c2)*{\theta}(t)=t5+tpto$

There is some error in my first edit and it should be the ${\int}_{- \infty}^{t}R_{55}(t-\tau){\theta}^{'}(\tau)d{\tau} $ And if when form $ -\infty $ to $0$ the ${\theta}(t) $ is undefined or is $0$ then It is ${\int}_{0}^{t}R_{55}(t-\tau){\theta}^{'}(\tau)d{\tau} $ And also by using the FindFit I get the similar r55,a more simple form so is there a way to solve ?

Remove["Global`*"];
omega = 1.2;

i1 = 25640;
ia = 4*10^6;
c2 = 1.4513*10^5;
t5 = 16000*Cos[omega*t];
tpto = 2.8*10^6*Theta'[t];

(*r55=1/(Pi (1.` +t^2)^3) 2 \
(7.200000000000001`*^7-2.1600000000000006`*^8 \
t^2+(-199396.49151683348` -317073.10946119425` \
t^2-130751.79771595637` t^4) Cos[10 t]+(238622.0308316204` \
t+388986.59820497024` t^3+163439.7471449455` t^5) Sin[10 t]);*)


r55[t_]=(44559034.3038137*t-26314556.6322484)/(-0.507157643287066-t^\
4);






s = NDSolve[{(i1 + ia)*Theta''[t] + 
 Integrate[r55*Theta'[Tau], {Tau, 0, t}] + c2*Theta[t] == 
t5 + tpto, Theta'[0] == 1, Theta[0] == 0}, Theta, {t, 0, 10}]
(*NDSolve::rdelay: Delay Tau is not real valued. >>*)
(*NDSolve::rdelay: Delay Tau is not real valued. >>*)



s = NDSolve[{(i1 + ia)*Theta''[t] + 
 NIntegrate[r55*Theta'[Tau], {Tau, 0, t}] + c2*Theta[t] == 
t5 + tpto, Theta'[0] == 1, Theta[0] == 0}, Theta, {t, 0, 10}]
(*NIntegrate::nlim: Tau = t is not a valid limit of integration. >>*)
(*NIntegrate::nlim: Tau = t is not a valid limit of integration. >>*)
(*NIntegrate::nlim: Tau = #1 is not a valid limit of integration. >>*)
(*General::stop: Further output of NIntegrate::nlim will be suppressed during this     calculation. >>*)
(*NDSolve::rdelay: Delay Tau is not real valued. >>*)
(*NDSolve::rdelay: Delay Tau is not real valued. >>*）

So is there a way that I can get the answer? Thanks!

Supplement

yes we can use the laplace transform to solve this question,bu the accuracy is the question,

  r55[t_]=(44559034.3038137*t-26314556.6322484)/(-0.507157643287066-t^\
4);

i get the laplace transform of r55,it is the function of p,so we call it r55p,but it is with Complex form，so i ge the fit of r55p,then i get

  r55p = 13829522.310484 + -3538064.23662678/p + 63545.307621876*p^2 + 
   235429.331970918/p^2 - 1546395.29415638*p;

and finally i get the solution of theta,it is

  Theta[t_] = 
  0.06541197518272371` E^(-33.8394608216844` t) - 
   0.38705195407627524` E^(-5.468703164047347` t) - 
   0.4851421230004468` E^(0.10884574410374667` t) + 
   0.8072987102672268` E^(
    0.18393223395811042` t) - (0.0005166083732284499` + 0.` I) Cos[
     1.2` t] + (0.000959271400198915` + 0.` I) Sin[1.2` t];


Plot[Theta[t], {t, 0, 10}]

and so we get the plot which is different from the answer 1

Answer 1

Completely revised Laplace-transform solution

In principle, this integro-differential equation can be solved by means of a Laplace transform, because it is linear and its integral is a convolution. However, the Laplace transform for r55 has a LeafCount of order 56000. Even the simplified version of r55, also presented in the question, is quite lengthy and involves derivatives of Hypergeometric functions. If, however, the Laplace transform and its inverse are performed numerically to obtain poles and residues, then considerable progress can be made. First, formally obtain the Laplace transform of the integro-differential equation.

Collect[LaplaceTransform[(i1 + ia)*Theta''[t] + 
    Integrate[r55[t - Tau]*Theta'[Tau], {Tau, 0, t}] + 
    c2*Theta[t] - (t5 + tpto), t, s], 
LaplaceTransform[Theta[t], t, s], Simplify] /. {Theta[0] -> 0, Theta'[0] -> 1};
sol = Solve[% == 0, LaplaceTransform[Theta[t], t, s]][[1, 1]] // Values
(* (4.02564*10^6 + (16000. s)/(1.44 + s^2))/
   (145130. - 2.8*10^6 s + 4.02564*10^6 s^2 + 1. s LaplaceTransform[r55[t], t, s]) *)

where LaplaceTransform[r55[t], t, s] is NIntegrate[r55 Exp[-s t], {t, 0, Infinity}] and is yet to be evaluated. To begin determining the inverse Laplace transform, first consider the poles represented by (16000. s)/(1.44 + s^2), which arise from the inhomogeneous term in the equation. Use Residue to obtain the residue of the numerator for each pole, divide them by the numerically evaluated denominators, and add them.

Residue[Numerator[sol] // Last, {s, 1.2 I}]/
    ReleaseHold[Denominator[sol] /. LaplaceTransform[r55[t], t, s] -> 
    Hold[NIntegrate[r55 Exp[-s t], {t, 0, Infinity}]] /. s -> 1.2 I]
term0 = ComplexExpand[2 Re[% Exp[I 1.2 t]]] // Chop
(* 0.000258352 - 0.000290531 I *)
(* 0.000516704 Cos[1.2 t] + 0.000581062 Sin[1.2 t] *)

Next, systematically obtain roots of the denominator numerically, again determine the corresponding residues, and for pairs of poles, sum them. (Residues are determined numerically by Cauchy's theorem, as described here.) The dominant term turns out to be

s /. FindRoot[145130.` - 2.8`*^6 s + 4.02564`*^6 s^2 + 
    1.` s NIntegrate[r55 Exp[-s t], {t, 0, 10^3}], {s, .05 + 0.05 I}, 
    Evaluated -> False]
Quiet@NIntegrate[
    sol /. LaplaceTransform[r55[t], t, s] -> NIntegrate[r55 Exp[-s t], {t, 0, 10^3}], 
    {s, % + eps (1 + I), % + eps (-1 + I), % + eps (-1 - I), 
     % + eps (1 - I), % + eps (1 + 1 I)}]/(2 Pi I)
term1 = ComplexExpand[2 Re[% Exp[%% t]]] // Chop
(* 0.0559559 + 0.0552868 I *)
(* 0.0992945 - 1.65608 I *)
(* 0.198589 E^(0.0559559 t) Cos[0.0552868 t] + 
   3.31215 E^(0.0559559 t) Sin[0.0552868 t] *)

The next largest oscillatory term is obtain in the same way.

(* 0.0405067 E^(-0.0140449 t) Cos[0.0059925 t] - 
   0.00359886 E^(-0.0140449 t) Sin[0.0059925 t] *)

In fact, there appears to be an infinite set of such slightly damped oscillatory terms, with frequencies near -.014 + 0.062 n I and their conjugates, with n an integer`. (I have computed ten of them.) There also is a weakly damped pure exponential.

s /. FindRoot[145130.` - 2.8`*^6 s + 4.02564`*^6 s^2 + 
    1.` s NIntegrate[r1 Exp[-s t], {t, 0, 10^3}], {s, -.014 + 0 I}, 
    Evaluated -> False];
Quiet@NIntegrate[
    sol /. LaplaceTransform[r[t], t, s] -> NIntegrate[r1 Exp[-s t], {t, 0, 10^3}], 
    {s, % + eps (1 + I), % + eps (-1 + I), % + eps (-1 - I), % + 
     eps (1 - I), % + eps (1 + 1 I)}]/(2 Pi I) // Chop;
 term13 = % Exp[%% t]
 (* 0.0201212 E^(-0.0141239 t) *)

Now, plot the sum of these terms.

Plot[{term0 + term1 + term2 + ... + term13}, 
    {t, 0, 10}, ImageSize -> Large, AxesLabel -> {t, Theta}, LabelStyle -> 
    Directive[Bold, Black, Medium]]

Comparing this plot with that in my other, purely numerical, answer shows a noticeable discrepancy for t < 3 but increasingly good agreement at larger t. A heavily damped oscillatory term appears to be missing here, probably with a frequency of order -1.3 + 4.5 I (and its conjugate). I have tried hard to find that term by the methods above but without success. Thus, the partly symbolic solution obtained here is valid only for larger t.

Answer 2

The problem proposed in the question above can be solved numerically by replacing the integro-differential equation by its finite difference counterpart, obtained in a straightforward fashion.

nm = 1000;
d = 10/nm;
th = ConstantArray[0, nm + 1]; th[[2]] = d;
dth = ConstantArray[0, nm]; dth[[1]] = d;
r = Reverse@Table[r55 /. t -> (n - .5) d, {n, 1, nm}];
Do[conv = dth[[;; n - 2]].r[[nm + 3 - n ;;]];
    th[[n]] = ((2.875457142857143` + 1.` d) th[[n - 2]])/(-2.875457142857143` + 1.` d) +
    ((5.750914285714287` - 0.10366428571428572` d^2) th[[n - 1]])/
        (2.875457142857143` - 1.` d) + 
    (7.142857142857143`*^-7 d^2 (conv - 16000.` Cos[1.2 (n - 1) d]))/
        (-2.875457142857143` + 1.` d);
    dth[[n - 1]] = th[[n]] - th[[n - 1]], {n, 3, nm + 1}]
plt = ListLinePlot[th, DataRange -> {0, 10}, PlotStyle -> Red, ImageSize -> Large, 
    AxesLabel -> {t, Theta},  LabelStyle -> Directive[Bold, Black, Medium]]

This solution may appear to be a simple exponential at late times. However, running it to, say t == 100, indicates that it is an exponentially growing oscillation at late times.

Also, the numerical procedure can be validated by solving the integro-differential equation with an expression for r55 that yields a symbolic solution by Laplace transform, such as 4.5 10^7 Cos[2.5 t] Exp[-1.1 t], which results in

-1.18557 E^(-0.388664 t) + 1.27708 E^(-0.04152 t) - 
    0.000880185 Cos[1.2 t] - 0.090631 E^(-0.537137 t) Cos[4.0469 t] + 
    0.00070688  Sin[1.2 t] + 0.134104 E^(-0.537137 t) Sin[4.0469 t]

and comparing it with the corresponding numerical solution. Results are identical to the eye. In general, these computations take only several seconds each.