Plot for pressure rise per wavelength

Question

I have another question that i need to ask. I have been struggling in making the graphs for pressure rise per wavelenth. my equation for pressure rise per wavelength is:

$Δp=∫_0^¹∫_0^¹(((∂p)/(∂x)))_{y=0}dt dx.$

and where

$(∂p)/(∂x)= (∂/(∂y))(s_{xy})-N²((∂ψ)/(∂y))+Grθ+Brσ$

and

$s_{xy}=ψ_{yy}[1-ζ(ψ_{yy})²]⁻¹. $

Now I have tried this code

F = Θ + a*Sin[2*π*(x - t)] + 
   b*Sin[2*π*(x - t) + ϕ];
N1 = Sqrt[M^2 + (1/k)];
x = 0.4;
t = 0.2;
a = 0.3;
b = 0.4;
m = 0.25;
ϕ = 2 π/3;
M = 1;
k = 0.8;
Nt = 0.8;
Nb = 0.4;
Pr = 0.2;
L = 0.1;
Bm = 4;
Bh = 2;
Θ = 1.5;
Gr = 0.7;
Br = 0;
ζ = 0.002;
Rn = 0.6;
h1 = -1 - m*x - a*Sin[2*π*(x - t) + ϕ]
h2 = 1 + m*x + b*Sin[2*π*(x - t)]
sol = NDSolve[{ψ''''[
      y] + ζ*(6 ψ''[y] *ψ'''[y]*ψ'''[y] + 
        3 ψ''[y]*ψ''[y]* ψ''''[y]) - 
     N1*N1*ψ''[y] + Gr*θ'[y] + Br*σ'[y] == 
    0, (1 + Pr*Rn)*θ''[y] + Nb*Pr*σ'[y]*θ'[y] + 
     Nt*Pr*θ'[y]*θ'[y] == 
    0, σ''[y] + Nt/Nb*θ''[y] == 0, ψ[h2] == F/
    2, ψ[h1] = -(F/2), ψ'[h1] == 0, ψ'[h2] == 
    0, σ'[h1] == Bm*σ[h1], σ'[h2] == 
    Bm*(1 - σ[h2]), θ'[h1] == 
    Bh*θ[h1], θ'[h2] == 
    Bh*(1 - θ[h2])}, {ψ, θ, σ}, {y, h1, h2}]
A1 = Plot[
  Evaluate[Integrate[((D[ψ[y], {y, 3}] + ζ*
          D[ψ[y], {y, 3}]*D[ψ[y], {y, 2}]*
          D[ψ[y], {y, 2}] + 
         2*D[ψ[y], {y, 2}]*D[ψ[y], {y, 2}]*
          D[ψ[y], {y, 2}]*D[ψ[y], {y, 3}] - 
         N1*N1*D[ψ[y], y] + Gr*θ[y] + Br*σ[y]) /. 
       y -> 0), {x, 0, 1}, {t, 0, 1}] /. sol], {y, h1, h2}, 
  PlotRange -> All, 
  PlotStyle -> {Darker[Blue, 0.5], Thickness[0.004]}, 
  AxesOrigin -> Automatic, 
  BaseStyle -> {FontFamily -> "Times", FontSize -> 15}, 
  FrameLabel -> {"Θ", "Δp"}, 
  Frame -> True, Axes -> False]

But it is giving me blank plot And errors are something like this :

NDSolve::deqn: "Equation or list of equations expected instead of -0.851076 in the first argument " , Integrate::ilim: Invalid integration variable or limit(s) in {0.4,0,1}. >> , ReplaceAll::reps: {NDSolve[{0.7 (θ^′)[y]-2.25 (ψ^′′)[y]+(ψ^(4))[y]+0.002 (Times[<<3>>]+Times[<<3>>])==0,0.16 (<<1>>^(<<1>>))[<<1>>]^2+0.08 (θ^′)[y] (σ^′)[y]+1.12 (θ^′′)[y]==0,2. (θ^′′)[y]+(σ^′′)[y]==0,ψ[1.48042]==0.851076,-0.851076,(ψ^′)[-1.03763]==0,(ψ^′)[1.48042]==0,(σ^′)[-1.03763]==4 σ[-1.03763],(σ^′)[1.48042]==4 (1-σ[<<1>>]),(θ^′)[-1.03763]==2 θ[-1.03763],(θ^′)[1.48042]==2 (1-θ[<<1>>])},{ψ,θ,σ},{y,-1.03763,1.48042}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >> NDSolve::dsvar: -1.03758 cannot be used as a variable. >>

ReplaceAll::reps: {NDSolve[<<1>>]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

NIntegrate::itraw: Raw object 0.4` cannot be used as an iterator. >>

NDSolve::dsvar: -1.03758 cannot be used as a variable. >>

and many more errors like these..Please help me out.

Answer 1

There are two problems in your code. First, there is a syntax error: ψ[h1] = -(F/2) instead of ψ[h1] == -(F/2) in the boundary conditions of the equation. After fixing this NDSolve returns a solution.

The second one is a wrong way you treat the expression under the integral. I take it you need to have this expression with the functions psi, sigma and theta taken from the solution. Then you need to make the substitution at once, before you put y->0, as below:

    ((D[ψ[y], {y, 3}] + ζ*D[ψ[y], {y, 3}]*
      D[ψ[y], {y, 2}]*D[ψ[y], {y, 2}] + 
     2*D[ψ[y], {y, 2}]*D[ψ[y], {y, 2}]*D[ψ[y], {y, 2}]*
      D[ψ[y], {y, 3}] - N1*N1*D[ψ[y], y] + Gr*θ[y] + 
     Br*σ[y]) /. sol /. y -> 0)

(*  {-2.74264}  *)

The result is a list with one element, a constant. I do not know the problem origin and, thus, cannot judge, if it is right, but I guess that it is not what you want. Further you integrate this number over x and t from 0 to 1, which is useless and gives the same constant. Let us note that from the very beginning nothing in your equation or in the expression under the integral depends upon x and t. Then you plot this constant from h1 to h2 which will give a straight line (of course, if you will take care of first removing the curly braces).

To summarize, the equation can be solved after removing a small syntax error, but the integral contains a conceptual mistake, and seems to be senseless in its present form.

Answer 2

The first problem is due to a typo. In the list of equations you feed to NDSolve you state

 ψ[h1] = -(F/2)

instead of

  ψ[h1] == -(F/2)

The next problem is the integration over x and t - first of all this is probably not what you want as you integrate over an object that does not depend on neither x nor t as you assigned x and t fixed values.

You might want to look at `ParametericNDSSolve', though, as this allows you to solve differential equations that still have an unevaluated parameter (e.g. and x and t).