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.

Improve this question
$\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
    Commented 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$
    – Sachin Kaushik
    Commented Dec 14, 2021 at 8:22
  • $\begingroup$ For your consideration, i have removed the table and updated the question. $\endgroup$
    – Sachin Kaushik
    Commented Dec 14, 2021 at 8:23
  • 2
    $\begingroup$ You give equations with machine precision (approx. 16) and then define WorkingPrecision->30. $\endgroup$
    – Daniel Huber
    Commented Dec 14, 2021 at 9:30

1 Answer 1

Reset to default
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   *)
Improve this answer
$\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$
    – Sachin Kaushik
    Commented 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$
    – Sachin Kaushik
    Commented 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
    Commented 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$
    – Sachin Kaushik
    Commented 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
    Commented Dec 15, 2021 at 12:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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