0
$\begingroup$

First I solved coupled ODE using NDSolve. Now I want to integrate the expression (p'''[0] - ((1.0/k)*p''[0]) + (beta*m*m*s[0]) with respect to the parameter x over the interval [0,1]. I need a plot of the integrated expression, say pr as a function of the parameter F. I am trying with the following codes.

a = 0.5;
b = 2.0;
t = 0.0;
Q0 = 1.0;
alp = 0.5;
alpha = 0.5;
beta=1.0;
k = 0.5;
m = 10.0;
h = 1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
eq1 = {s''[y] == (m/alpha)*Sinh[alpha*s[y]]};
bc1 = {s'[0] == 0, s[h] == 1};
eq2 = {p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0};
bc2 = {p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0};
px = ParametricNDSolve[
   Join[eq1, eq2, bc1, bc2], {s, p}, {y, 0, h}, {x, F} ];
pr = Evaluate[
  NIntegrate[(p'''[y] - ((1.0/k)*p''[y]) + (beta*m*m*s[y])), {x, 0, 
     1}] /. {y -> 0}]

The graph of Integrated value vs F should look like the following graph. Kindly note that here beta1 is another parameter involved in the problem. But for simplicity, I ignored that parameter. Here we can try for any other parameter that is present in the code itself, for example, the parameter m. Graph

$\endgroup$
6
  • $\begingroup$ There are several similar questions posted by you. Instead of giving feedback you prefer to post new questions, no good style! $\endgroup$ Sep 14, 2023 at 8:44
  • $\begingroup$ I am not getting the desired answer. That's why I am posting. If you are unable to help then please ignore. Thank you. $\endgroup$ Sep 14, 2023 at 8:55
  • $\begingroup$ If you "don't get the desired answer" you should give a feedback. This interaction makes stackexchange work and honors the effort made by users trying to help you! $\endgroup$ Sep 14, 2023 at 9:04
  • $\begingroup$ O.K. Thank you. $\endgroup$ Sep 14, 2023 at 9:10
  • 1
    $\begingroup$ How is parameter \[Beta] defined? $\endgroup$ Sep 14, 2023 at 9:32

1 Answer 1

1
$\begingroup$

Your definitions (keep h undefined)

ClearAll["Global`*"]
a = 0.5;
b = 2.0;
t = 0.0;
Q0 = 1.0;
alp = 0.5;
alpha = 0.5;
beta = 1.0;
k = 0.5;
m = 10.0;

eq1 = {s''[y] == (m/alpha)*Sinh[alpha*s[y]]};
bc1 = {s'[0] == 0, s[h] == 1};
eq2 = {p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0};
bc2 = {p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0};

function ph[...] returns p'''[y] - (1/k)*p'[y] + (beta*m*m*s[y])

ph[h_, F_] :=Block[{sol}, 
sol = NDSolveValue[{s''[y] == (m/alpha)*Sinh[alpha*s[y]],p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0, 
s'[0] == 0,s[h] == 1, p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0}, 
p'''[y] - (1/k)*p'[y] + (beta*m*m*s[y]), {y, 0, 1} ]    
]

integration

xi = Range[0, 1, .02];
hi = Map[1.0 - (a*Cos[Pi*(# - t)]*Cos[Pi*(# - t)]) &, xi];
data = MapThread[{#1, NIntegrate[ph[#2, 1], {y, 0, 1}]} &, {xi, hi}];
ListPlot[data]

enter image description here

addendum( modified with new function ph0[] )

ph0[h_, F_] := 
Block[{sol}, 
sol = NDSolveValue[{s''[y] == (m/alpha)*Sinh[alpha*s[y]],p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0, s'[0] == 0,s[h] == 1, p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0}, 
p'''[0] - (1/k)*p'[0] + (beta*m*m*s[0]), {y, 0, 1}]]


Fi = Range[-3, 3];
Fdp=Map[{#, NIntegrate[ph0[1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]), #],{x, 0,1}]} &, Fi    ]//Quiet
ListPlot[{Fdp, Fdp}, AxesLabel -> {F, dp}, Joined -> {True, False}]

enter image description here

$\endgroup$
15
  • $\begingroup$ It seems to be fine. But here you have considered the range of the dependent variable y as [0,1]. Also as I mentioned in my question I want to plot integration of (p'''[0] - ((1.0/k)*p''[0]) + (beta*m*m*s[0]) with respect to the variable x vs the parameter F. However, the desired graph should look different. I have edited the question. Please see. $\endgroup$ Sep 14, 2023 at 9:22
  • $\begingroup$ Thank you very much. This is the required result. But when I try to run it in my system the code is running for 3-4 hours with lots of errors and I don't get any output. $\endgroup$ Sep 15, 2023 at 10:31
  • $\begingroup$ See my modified answer. Sorry, don't know we it evaluated yesterday $\endgroup$ Sep 15, 2023 at 11:05
  • $\begingroup$ Thank you. Now it is working fine. However, I've two questions for my understanding. 1. Why do you consider the range of y, [0,1] instead of [0,h]? 2. What is the reason for writing {Fdp,Fdp} twice in ListPlot[{Fdp, Fdp}, AxesLabel -> {F, dp}, Joined -> {True, False}]? $\endgroup$ Sep 15, 2023 at 11:50
  • $\begingroup$ Also, is it possible to plot the variation for four different values of the variable m say m=5, 10, 15, 20 in the same figure? $\endgroup$ Sep 15, 2023 at 12:05

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.