@Fateh
You asked (1)
explain the purpose of the command
rhs(sol[j](0.5)[3])
In the code I supplied, I solve your ode for different entries of values n the list 'params' by looping through the entries in the list [0.25, 0.3, 0.4, 0.5, 0.7], for values of 'N', while holding 'M' and 'K' fixed at 0.2 and 2 respectively. The solution for each value of the loop index 'j' is assigned to sol[j] - so for example sol[4] is the solution corresponding to the parameter list [M=0.2, N=0.5, k=2].
What the dsolve() command returns is actually a procedure, to which one can supply an argument and obtain values - so for example sol[4](0.5) will return a list of values for your system when x=0.5. If you try this in my worksheet, it will return
[x = 0.5, p(x) = 0.108665365029661, lambda = 0.658681811222779]
In other words the value of 'x', the value of 'p(x)' (with x=0.5) and the value of lambda. sol[4](0.5)[3] will return the third entry of this list, which is lambda = 0.658681811222779; rhs(sol[4](0.5)[3]) just takes the right-hand-side of the supplied expression so will return 0.658681811222779. Thus
rhs(sol[j](0.5)[3])
is just a "short" way to get the value of lambda from the j-th solution, for the purposes of printing. Note that there is no particular reason for uing the value x= 0.5 in the above - any value between 0.0 and 1.0 would give the same answer, since lambda is constant for each solution. The value of lambda varies (slightly) when one changes the parameter set, so the command
seq( printf( "%12.8f %12.8f\n", par[j], rhs(sol[j](0.5)[3])), j=1..numelems(par));
prints out the j-th value of list [0.25, 0.3, 0.4, 0.5, 0.7], which is being used for 'N' and the corresponding value of lambda from the j-th solution
You asked (2)
How could i change my method from BVP to rk45?
You can't! Runge-Kutta methods such as rk45, can only be used on initial value problems; you have a boundary value problem for which differeent solution methods are necessary. It is generally possible to convert a boundary value problem to an initial value problem using something called the "shooting method", but there seems little point adding this extra complexity when Maple's BVP solution methods are working well
You asked (3)
I want to solve the given first order ode with mentioned boundary condition for the expression of p(x) then I need to plot p(x) for the different involved parameters.
This is exactly what my previous worksheet does, it uses the parameter lists
[M=0.2, N=0.25, k=2]
[M=0.2, N=0.3, k=2]
[M=0.2, N=0.4, k=2]
[M=0.2, N=0.5, k=2]
[M=0.2, N=0.7, k=2]
plots the curve p(x) for each of the above parameter lists, and prints the value of lambda. Obviously you can use any lists of parameters which you want. The above was just an arbitrary choice on my part
You asked (4)
I want to calculate the expression and table for lambda for different values of involved parameters by fixing some and changes one.
Since everyting is being calculated numerically you will not be able to obtain an "expression" for lambda. In the worksheet I supplied previously, you did get a "table" of lambda . So you will always be able to calculate a value for lambda for any supplied parameter list, ie any given value for M, N, k
You asked (5)
next step is to find the maximum value of p(x) and the point where it is maximum then the plot for pmax?
This one is interesting! It is actually quite trivial to obtain the maximum value of p(x) and then produce plots of p(x)/max(p(x)), whic is what *appears* to be plotted in Figure 2 of the pdf whixh you supplied. When I did this for the parameter sets listed above - the curves came out as identical. In other words the parameter N is basically a "scale factor" whose effect is completely cancelled when one rescales the curves using the appropriate max(p(x)) value.
For this reason in the attached new worksheet,I fixed the parameters 'M', 'N' and varied 'k' instead, so it produces outputs for the combinations
[M=0.2, N=0.4, k=1.0]
[M=0.2, N=0.4, k=2.0]
[M=0.2, N=0.4, k=3.0]
[M=0.2, N=0.4, k=4.0]
[M=0.2, N=0.4, k=5.0]
This now produces different plots of p(x)/max(p(x)) for each parameter set - and they look (roughly!) like the curves in Figure 2 of the supplied pdf
You asked (6)
At the last again want to integrate p(x) from 0 to 1 and generate a table for load for different values of the involved parameter.
Not too difficult to produce these integrals, although I'm uncertain whether you want to produce the integral of p(x) or the integral of p(x)/max(p(x)), so the following code and the attached worksheet odeProb3.mw does both
restart;
h:= k - (k - 1)*x;
DE:= diff(p(x), x)
= -2*M^2
*
( -exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
+ M*sinh(M*h/sqrt(N))*lambda + exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)
)
/
( exp(-M*h/sqrt(N))
*
( M*exp(2*M*h/sqrt(N))*h
- 2*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
- M*h + 4*exp(M*h/sqrt(N))*sqrt(N)
- sqrt(N)*exp(2*M*h/sqrt(N))
- sqrt(N)
)
);
BC := p(0) = 0, p(1) = 0;
# par:=[0.25, 0.3, 0.4, 0.5, 0.7];
par:= [ 2.0, 3.0, 4.0, 5.0, 6.0];
colors:=[red, blue, green, cyan, black]:
for j from 1 by 1 to numelems(par) do
params := [M = 0.2, N = 0.4, k =par[j]];
#
# Solve the ode for parameter values in the
# list called params above
#
sol[j] := dsolve
( eval
( {BC, DE},
params
),
numeric,
abserr=1e-03
);
#
# Compute the x-value at which the maximum
# value of p(x) occurs
#
xmax[j]:= fsolve
( evalf
( rhs
( eval
( DE,
[ [M = 0.2, N = 0.4, k =par[j]][], sol[j](0.5)[3]]
)
)
),
x=0..1
);
#
# At the obtained x-value above, determine what
# the maximum value of p(x) actually is
#
pmax[j]:= rhs(sol[j](xmax[j])[2]);
#
# Plot p(x) versus x
#
plt1[j]:= plots:-odeplot
( sol[j],
[x, p(x)],
x = 0 .. 1,
color=colors[j]
);
#
# Plot p(x)/max(p(x)) versus x
#
plt2[j]:= plots:-odeplot
( sol[j],
[x, p(x)/pmax[j] ],
x = 0 .. 1,
color=colors[j]
);
od:
#
# Display all plots of p(x) versus x with some
# annotations
#
plots:-display( seq
( plt1[j],
j=1..numelems(par)
),
title=typeset( p(x), "versus ", x),
titlefont=[times, bold, 16],
legend=[ seq
( typeset("k= ", par[j]),
j=1..numelems(par)
)
]
);
#
# Display all plots of p(x)/max(p(x)) versus x with some
# annotations
#
plots:-display( seq
( plt2[j],
j=1..numelems(par)
),
title=typeset( p(x)/pmax, "versus ", x),
titlefont=[times, bold, 16],
legend=[ seq
( typeset("k= ", par[j]),
j=1..numelems(par)
)
]
);
#
# Print an output table of lambda as k varies
#
printf( "\n\n%12s %12.8s\n", "k", "lambda");
seq( printf( "%12.8f %12.8f\n", par[j], rhs(sol[j](0.5)[3])), j=1..numelems(par));
#
# Print an output table of the integral of p(x) and
# the integral of P(x)/max(p(x)) as k varies
#
with( Student[Calculus1] ):
printf( "\n\n%8s %12s %20s\n",
"k",
"Int( p(x) )",
"Int( p(x)/max(p(x)) )"
);
f1:= proc(x, j)
if type(x, numeric)
then return rhs(sol[j](x)[2]);
else return 'procname(_passed)';
fi;
end proc:
f2:= proc(x, j)
if type(x, numeric)
then return rhs(sol[j](x)[2])/pmax[j];
else return 'procname(_passed)';
fi:
end proc:
seq( printf( "%12.8f %12.8f %12.8f\n",
par[j],
ApproximateInt( f1(x, j), x=0..1, partition=99, output=value),
ApproximateInt( f2(x, j), x=0..1, partition=99, output=value)
),
j=1..numelems(par)
);
This produces the graphs
and the tables
k lambda
2.00000000 0.65672680
3.00000000 0.73526777
4.00000000 0.78216130
5.00000000 0.81286188
6.00000000 0.83448106
k Int( p(x) ) Int( p(x)/max(p(x)) )
2.00000000 0.06184010 0.63537105
3.00000000 0.05711187 0.59111658
4.00000000 0.04762111 0.55113440
5.00000000 0.03954397 0.51673786
6.00000000 0.03322128 0.48724002
