# How to solve this differential equation with integral?

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

• Theta[Tau] is undefined for Tau < 0, and NIntegrate should be used instead of Integrate. However, even with these items corrected, it is not clear that NDSolve can solve such an integro-differential equation. By the way, I suggest you remove the first two blocks of code from the question, because they do nothing. I also suggest that you include the first few error messages that you receive. – bbgodfrey Mar 19 '18 at 2:09
• Incidentally, DSolve with Integrate can, in principle, solve such integro-differential equations, but not this one. – bbgodfrey Mar 19 '18 at 2:23
• @bbgodfrey r55 is definded in braket in the code of first part. – dcydhb Mar 19 '18 at 13:50
• Actually, not. The code contains two comments, each including a possible definition of r55, but as comments neither is available to the code that follows. Note also that the comments contain r55[t] and the subsequent lines of code just r55, which is not the same. To obtain useful advice from readers, please display the code you actually used and what error messages, if any, you actually received. In the event of many consecutive error messages, including just the first few typically is sufficient. – bbgodfrey Mar 19 '18 at 14:05
• It seems R55 is not actually a function of tau, so you can factor out of the integral.. (Or should the t in the R55 expression really be t-tau ? ) – george2079 Mar 19 '18 at 14:40

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.

• ,thanks for your solution!!! but the difference of the plot of r1,r2,r3when t=2may make a great difference. – dcydhb Mar 25 '18 at 3:46
• Using the Riemann Merlin inversion formula is a novel method, but how to find all the poles is crucial, thank you very much! – dcydhb Apr 5 '18 at 8:19
• I believe that I have found all the poles with significant contributions to the solution, except the one that I mentioned in the answer. Additionally, there may be branch cuts, although none were obvious, when I plotted parts of the complex plane. Incidentally, the integral from zero to infinity of your original r55 is zero, your approximation to it about - 500000, and of my original approximation about 7000000, which led to excessive damping of my earlier solution, as you correctly identified. – bbgodfrey Apr 5 '18 at 13:01
• Incidentally, the dominant asymptotic term of the solution for your simplified approximation to r55 is 0.235231 E^(0.0576946 t) Cos[0.0467374 t] + 3.59494 E^(0.0576946 t) Sin[0.0467374 t]. – bbgodfrey Apr 5 '18 at 13:09
• Thank you for your efforts, I think you have solved my problem very effectively, thank you very much! – dcydhb Apr 6 '18 at 8:42

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.