2
$\begingroup$

I'm trying to get solution of the following ode using following code. But no plotting appears. I need help to solve the problem.

a = 0.5;
b = 2.0;
x = 0.1;
t = 0.1;
Ha = 1.0;
Q0 = 1.0;
we = 0.2;
alp = 1.0;
F = Q0*Exp[-alp*t]
dVals = {0.01, 0.1, 0.2, 0.3};
h = 1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
sys2 = {(1.0 - d) s''''[y] + (2.0*d*we*s''[y]*s''''[y]) + (2.0*d*we*
       s''[y]*s''[y]) == Ha*Ha*s'[y]};
iv2 = {s[-h] == -F/2.0, s[h] == F/2.0, s'[-h] == 0.0, s'[h] == 0.0};
sol2[d_] := 
  NDSolve[Join[sys2, iv2], s, {y, -h, h}, 
   Method -> {"StiffnessSwitching", 
     Method -> {"ExplicitRungekutta", Automatic}}, 
   WorkingPrecision -> 35, AccuracyGoal -> 5, PrecisionGoal -> 4];
Plot[Evaluate[Table[s'[y] /. sol2[d], {d, dVals}]], {y, -1, 1}, 
 PlotRange -> Full, AxesLabel -> {y, \[CapitalPsi]}, Frame -> False, 
 AxesStyle -> Black, 
 PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted, 
    Black}, {Dashed, Black}}, 
 PlotLegends -> 
  Placed[{"d = 0.01  ", "d = 0.1  ", "d = 0.2  ", "d = 0.3 "}, 
   Center]]
$\endgroup$
2
  • $\begingroup$ we is undefined. $\endgroup$
    – bbgodfrey
    Jul 19, 2017 at 12:47
  • $\begingroup$ Value of we is very small. It may be take as 0.2. I edited the question. Hope it will be okay to get a solution now. $\endgroup$ Jul 19, 2017 at 13:19

1 Answer 1

3
$\begingroup$

It is good practice to use rational rather than decimal numbers with NDSolve when possible. Also, the equation actually is not stiff, despite the error message that sometimes occurs. So, use

a = 1/2; b = 2; x = 1/10; t = 1/10; Ha = 1; Q0 = 1; alp = 1; F = Q0*Exp[-alp*t]; we = 2/10;
dVals = {1/100, 1/10, 2/10, 3/10};
h = 1 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
sol2[d_] := Flatten@NDSolve[{(1 - d) s''''[y] + (2*d*we*s''[y]*s''''[y]) + (2*d*
   we*s''[y]*s''[y]) == Ha*Ha*s'[y], s[-h] == -F/2, s[h] == F/2, s'[-h] == 0, s'[h] == 0}, 
   {s, s', s'', s'''}, {y, -h, h}, WorkingPrecision -> 30]
Plot[Evaluate[Table[s'[y] /. sol2[d], {d, dVals}]], {y, -h, h}, 
    PlotRange -> Full, AxesLabel -> {y, Ψ}, Frame -> False, AxesStyle -> Black, 
    PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted, Black}, {Dashed, Black}}, 
    PlotLegends -> Placed[{"d = 0.01 ", "d = 0.1 ", "d = 0.2 ", "d = 0.3 "}, Center]]

enter image description here

Note, also, that the Plot range has been reduced to correspond to the NDSolve range.

Addendum

It may seem strange that s'[t] is essentially unchanged as d is increased from near zero to 3/10. In fact, s[t] and s''[t] also are essentially unchanged. s'''[t], on the other hand changes significantly, as shown here.

Plot[Evaluate[Table[s'''[y] /. sol2[d], {d, dVals}]], {y, -h, h}, 
    PlotRange -> Full, AxesLabel -> {y, Ψ}, Frame -> False, AxesStyle -> Black, 
    PlotStyle -> {{Line, Black}, {DotDashed, Black}, {Dotted, Black}, {Dashed, Black}}, 
    PlotLegends -> Placed[{"d = 0.01 ", "d = 0.1 ", "d = 0.2 ", "d = 0.3 "}, Center]]

enter image description here

The third derivative is singular for y near h for d == 3/10, because the coefficient of s''''[y] vanishes there,

Coefficient[(1 - d) s''''[y] + (2*d*we*s''[y]*s''''[y]) + (2*d*we*
    s''[y]*s''[y]) - Ha*Ha*s'[y], s''''[y]]
(* 1 - d + 2/5 d s''[y] *)

Indeed, the coefficient vanishes for d == 3/10 and s''[h] approximately equal to -5.4. Empirically, the solution is on the verge of becoming singular at y == h for d == 2985209/10000000.

Plot[Evaluate[{s[y], s'[y], s''[y], s'''[y] + 10} /.   sol2[2985209/10000000]],
   {y, -h, h}, PlotLegends -> Placed[{s[y], s'[y], s''[y], s'''[y] + 10}, {.8, .8}]]

enter image description here

$\endgroup$
4
  • $\begingroup$ Thank you. I have to show how the solution is affected by the two key parameters 'd' and 'we' in two different plots. In this case no charges take place in the above plot. Then what is the advantages of taking such parameter in ode. I'm really confused. $\endgroup$ Jul 19, 2017 at 13:38
  • $\begingroup$ @BiswajitMallick With larger values of we, the answer is likely to change, although solving the ODE will be more difficult. $\endgroup$
    – bbgodfrey
    Jul 19, 2017 at 13:54
  • $\begingroup$ @BiswajitMallick I'll look more at this. $\endgroup$
    – bbgodfrey
    Jul 19, 2017 at 15:05
  • $\begingroup$ Thank you very much. @bbgodfrey $\endgroup$ Jul 19, 2017 at 16:03

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.