0
$\begingroup$
Table[sbar = 100;
theta = 1;
h = 5/10;
w2 = i;
 nu = 7/10;
al[t_] = sbar Sin[w2 t];
rho[t_] = 1 + nu Sin[w2 t];
r[x_, t_] = -theta (1 - x) D[rho[t], t] // Simplify;
 r1[x_, t_] = -theta x (1 - x) al[t] rho[t] (1 - 2 x + h (4 x - 3)) //
Simplify;

 (eqs = Subtract @@@ {D[v[x, t], 
        t] == -al[t] x (1 - x) D[(x + h (1 - 2 x)) v[x, t], 
          x] + (x (1 - x)/(2 rho[t])) D[v[x, t], x, x] + r[x, t] +
         r1[x, t], v[x, 0] == 0, v[0, t] == 0, v[1, t] == 0} // 
    Together // Numerator // Expand // Simplify) // TableForm;

 vsol = v /. 
  First@NDSolve[Thread[eqs == 0], v, {x, 0, 1}, {t, 0, 100 Pi/w2}(*,
 MaxSteps\[Rule]10^5,MaxStepSize\[Rule]{.0002,.04},
 StartingStepSize\[Rule]{.0001,.01}*), WorkingPrecision -> 30];
 nint1 = NIntegrate[2 vsol[x, t], {x, 0, 1}, {t, 52 Pi/w2, 54 Pi/w2}, 
MaxRecursion -> 20];
 Integrate[
2 (1 - x) rho[t], {x, 0, 1}, {t, 52 Pi/w2, 54 Pi/w2}]/(2 Pi/w2) + 
nint1/(2 Pi/w2), {i, 10, 90, 40}]

Please suggest me how to remove the numerical errors, clearly the answer which is coming is clearly off and absurd.

$\endgroup$
4
  • 3
    $\begingroup$ It would be helpful if you added a few more words of what you think the problem is and what you expected. Also, it might be useful for you to simplify your code. When I look at it and see a Table[...,{i,420,420,420}] I am not sure I want to spend time on this. $\endgroup$
    – user21
    Dec 14, 2021 at 7:26
  • $\begingroup$ Sir the answer was expected to be around 2. I wrote table like this because i needed values for various i. But for i larger than 400, my answers are not correct. You can ignore this table thing because i just wanted to calculate for one particular i. $\endgroup$ Dec 14, 2021 at 8:22
  • $\begingroup$ For your consideration, i have removed the table and updated the question. $\endgroup$ Dec 14, 2021 at 8:23
  • 2
    $\begingroup$ You give equations with machine precision (approx. 16) and then define WorkingPrecision->30. $\endgroup$ Dec 14, 2021 at 9:30

1 Answer 1

2
$\begingroup$

Take the periodicity of the diffequation into account. Since t only occurs in Sin or Cos with minimal 10 t, you have peroidicity of 2 Pi/10. NDSolve for high t would be very unreliable; integrate within lower period (But not from the beginning t==0 because of disturbing effects due to initial conditions). Rationalize parameters.

sbar = 420;
theta = 1;
h = 1/2;
w2 = 10;
nu = 7/10;
al[t_] = sbar Sin[2 w2 t];
rho[t_] = 1 + nu Sin[w2 t]; 
r[x_, t_] = -theta (1 - x) D[rho[t], t] // Simplify;
r1[x_, t_] = -theta x (1 - x) al[t] rho[t] (1 - 2 x + h (4 x - 3)) // 
 Simplify;

(eqs = Subtract @@@ {D[v[x, t], 
       t] == -al[t] x (1 - x) D[(x + h (1 - 2 x)) v[x, t], 
         x] + (x (1 - x)/(2 rho[t])) D[v[x, t], x, x] + r[x, t] + 
       r1[x, t], v[x, 0] == 0, v[0, t] == 0, v[1, t] == 0} // 
   Together // Numerator // Expand // Simplify) // TableForm

vsol = v /. 
  First@NDSolve[Thread[eqs == 0], v, {x, 0, 1}, {t, 0, 3}, 
MaxSteps -> 10^5, MaxStepSize -> {.0002, .04}, 
StartingStepSize -> {.0001, .01}]

Periodiity in t for all x shown in graphic. Integrating from t=2 Pi/w2 to t=4 Pi/w2 should yield the same as your given range.

Manipulate[
  Plot[{vsol[x, t], vsol[x, t + 2 Pi/w2], vsol[x, t + 4   Pi/w2]}, {t, 
  2 Pi/w2, 4 Pi/w2}, PlotPoints -> 100], {{x, 1/2}, 0, 1}]

Plot3D[vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}, PlotRange -> All,
 PlotPoints -> 50]

nint1 = NIntegrate[2 vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}]

(*   -0.00674974   *)

(*   Error message: ....NIntegrate obtained -0.00674974 and 6.36785*10^-8 for the integral \
and error estimates   *)

If you are not satisfied with the estimated error of 6 10^-8, you have to use higher WorkingPrecision. But this takes very long and may crash due to not enough memory.

Integrate[
 2 (1 - x) rho[t], {x, 0, 1}, {t, 42 Pi/w2, 44 Pi/w2}]/(2 Pi/w2) + 
nint1/(2 Pi/w2)

(*   0.989257   *)
$\endgroup$
5
  • $\begingroup$ Thank you so much!!!!!! but the periodicity for small t is not strong. I wanted to avoid the initial disturbances and thats why integrating over large time. I will look at it. Many thanks. $\endgroup$ Dec 15, 2021 at 6:59
  • $\begingroup$ Could you also please mention that while I am integrating over x and try to plot it for time only. It is not able to integrate . Plot[{NIntegrate[vsol[x, t], {x, 0, 1}]}, {t, 2 Pi/w2, 4 Pi/w2}] $\endgroup$ Dec 15, 2021 at 7:05
  • 1
    $\begingroup$ Tell NIntegrate not to try symbolicProcessing for that interpolating function vsol Plot[NIntegrate[vsol[x, t], {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}], {t, 2 Pi/w2, 4 Pi/w2}] $\endgroup$
    – Akku14
    Dec 15, 2021 at 7:57
  • $\begingroup$ Sir, could you please look the edited one, this also contains error and give me wrong solution. $\endgroup$ Dec 15, 2021 at 11:53
  • $\begingroup$ On my older version 8.0 this works well. As alternative you can try wsol[x_?NumericQ, t_?NumericQ] = vsol[x, t]; nint[t_?NumericQ] := NIntegrate[wsol[x, t], {x, 0, 1}]; Plot[nint[t], {t, 2 Pi/w2, 4 Pi/w2}] $\endgroup$
    – Akku14
    Dec 15, 2021 at 12:10

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.