## 20214 Reputation

15 years, 318 days

## Thank you...

@Carl Love  for the useful information on the binary search.

## ?...

@Carl Love  I did not understand the meaning of your comment. Did you mean that the procedure should check that the given list already is sorted and produce an error message otherwise? Or something else?

## Re...

You can leave the arrows in your code, but then you should write  g:=x->f(x)+a[i]*x^(k-1)  instead of  f:=x->f(x)+a[i]*x^(k-1) , because a contradiction arises: f  is determined through f. In fact, on the left you have one function, and on the right - the other.

This works:

restart;
P:= 7919; #(this is the1000-th prime)
S:=5000; k:=5:
a:=[0\$i=1..k];
f:=[0\$i=1..k];
a[1]:=S;
f:=x->a[1];

with(RandomTools):
for i from 2 to k do
a[i]:=Generate(integer(range=1..P-1));
g:=x->f(x)+a[i]*x^(k-1);
g(x);
od;

## textplot...

@ravenHound  May be  plots:-textplot  command will be usefull:

restart;
f:=(x,c,k)->c*x/(k+x):
plots:-display(plot([f(x,5,2), f(x,5,3)], x=0..50, color=[red,blue], size=[800,400]), plots:-textplot([[12,4.5,"k=2"], [12,3.8,"k=3"]], font=[times,16])) ;

## No errors...

I just copied your code into the Maple 2018 worksheet. There are no errors:

restart;
g := x -> (10*x)/(3+x);

f := x -> (10*x)/(5+x);

h := x -> (8*x)/(5+x);

plots:-display(
plot( [g(x),f(x), h(x)], x=0..30, color=[red,blue, green], legend = ["k=3,c=10", "k=5, c=10","k=5,c=8"]));

plot( [g(x),f(x), h(x)], x=0..30, color=[red,blue, green], legend = ["k=3,c=10", "k=5, c=10","k=5,c=8"]);

works exactly the same.

## Initial approximation...

@zdamasceno  Now everything is clear. At first I just did not understand the meaning of what you first called the boundary condition. In fact, this is just an initial approximation using Newton's method. In Maple, you do not need to specify it. Here are 2 solutions to your example. In the first, this initial approximation is indicated. In the second solution, the range for the wanted solution is specified instead. The results (lists  S1  and  S2) are identical.

 > restart; eqns := {(1-sin(b_u)*sin(s)/(cos(b_u-y_f)*cos(s-y_f)))*(1/cos(s-y_f)^2-1/sin(s)^2) = -(cot(s)+tan(s-y_f)+Z)*sin(b_u)*(cos(s)/cos(s-y_f)+sin(s)*sin(s-y_f)/cos(s-y_f)^2)/cos(b_u-y_f)}; b_u := 1/tan(0.8); Z := 892/(27417000*f_z); y_f := 9*Pi/180; b:=(1/4)*Pi-1/2*(b_u-y_f); S1:=[seq(rhs(fsolve(eval(eqns, f_z=a), s=b)[]), a=0.00005..0.0005, 0.000001)]; S2:=[seq(rhs(fsolve(eval(eqns, f_z=a), s=0..1)[]), a=0.00005..0.0005, 0.000001)]; plots:-implicitplot(eqns[], f_z=0..0.0005, s=0..0.5, gridrefine=5);
 >

## The meaning of the condition s_0 := (1/4...

@zdamasceno  For any value of the parameter  f_z  , you get a trigonometric equation for  s - variable, which has the infinite number of solutions. The condition   s= (1/4)*Pi-1/2*(b_u-y_f)  when f_z=0  allows us to choose some branch of solutions. But if we substitute  f_z=0  into the equation, we get an error (the divizion by 0). I guess that you have incorrectly set  Zf_z - variable must be in the numerator (not in the denominator). In this case we have  s= (1/4)*Pi-1/2*(b_u-y_f)=0.3783306793   for  f_z=0 :

restart;
eqns := {(1-sin(b_u)*sin(s)/(cos(b_u-y_f)*cos(s-y_f)))*(1/cos(s-y_f)^2-1/sin(s)^2) = -(cot(s)+tan(s-y_f)+Z)*sin(b_u)*(cos(s)/cos(s-y_f)+sin(s)*sin(s-y_f)/cos(s-y_f)^2)/cos(b_u-y_f)};
b_u := 1/tan(0.8);
Z := 892/27417000*f_z;
y_f := 9*Pi/180;
b:=fsolve(eval(eqns, f_z=0), s=0..1)[];
(1/4)*Pi-1/2*(b_u-y_f);
# This equals b
S:=[seq(rhs(fsolve(eval(eqns, f_z=a), s=0..1)[]), a=0.00005..0.0005, 0.000001)];

@zdamasceno  If I understood correctly, then just add 1 point  [0, (1/4)*Pi-1/2*(b_u-y_f)]:

restart;
eqns := {(1-sin(b_u)*sin(s)/(cos(b_u-y_f)*cos(s-y_f)))*(1/cos(s-y_f)^2-1/sin(s)^2) = -(cot(s)+tan(s-y_f)+Z)*sin(b_u)*(cos(s)/cos(s-y_f)+sin(s)*sin(s-y_f)/cos(s-y_f)^2)/cos(b_u-y_f)};
b_u := 1/tan(0.8);
Z := 892/(27417000*f_z);
y_f := 9*Pi/180;
S:=[(1/4)*Pi-1/2*(b_u-y_f), seq(rhs(fsolve(eval(eqns, f_z=a), s)[]), a=0.00005..0.0005, 0.000001)];
plot([0,seq(a, a=0.00005..0.0005, 0.000001)], S);

@Muhammad Usman  The following code produces the desired result:

```restart;

alpha := 1: k := 2: M := 3:
printlevel := 3;

for n from 1 while n <= 2^(k-1) do

for m from 0 while m <= M-1 do

for j from 0 while j <= M-1 do

Omega[m, j] := 2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+alpha)*GAMMA(alpha)^2/(Pi*2^(1-2*alpha)*GAMMA(j+2*alpha)))*(sum((-1)^i*GAMMA(j-i+alpha)*2^(j-2*i)*(sum((1/2)*binomial(m, l)*(2*n-1)^(m-l)*(1+(-1)^(j-2*i+l))*GAMMA((1/2)*j-i+(1/2)*l+1/2)*GAMMA(alpha+1/2)/GAMMA(alpha+1+(1/2)*j-i+(1/2)*l), l = 0 .. m))/(GAMMA(alpha)*factorial(i)*factorial(j-2*i)), i = 0 .. floor((1/2)*j)))/2^(k*(m+1));
A[n]:=Matrix(3, (i,j)->Omega[i-1,j-1]/sqrt(2)/sqrt(Pi));
od;  od;  od;
A[1];
A[2];
<A[1], Matrix(3); Matrix(3), A[2]>;
```

Edit.

## OK...

@waseem  Probably you have not read about what I wrote above: "Unfortunately, the developers did not provide an image of the points in the legends. This can be done manually, but rather cumbersome."

## Workaround...

@waseem

Probably you have an older version of Maple, which does not have this option. Try this code:

```restart;
h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=(4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4:
lambda:=(F,Nb,delta2)->Int(K1,z=0..1):

L1:=0.2:
d1:=0.2:
alpha:=Pi/6:
A:=plot( [seq(seq(lambda(F,Nb,delta2), Nb=[0.1,0.2,0.3]), F=[0.1,0.2,0.5])], delta2=0.02..0.1, style=line, linestyle = [solid,longdash,dashdot],     thickness = 2,color=[red\$3,blue\$3,black\$3], legend=["Nb=0.1",""\$2,"Nb=0.2",""\$2,"Nb=0.3",""\$2]):

B:=plot( [seq(seq(lambda(F,Nb,delta2), Nb=[0.1,0.2,0.3]), F=[0.1,0.2,0.5])], delta2=0.02..0.1, style=point, color=[red\$3,blue\$3,black\$3],       symbol=[solidbox, solidcircle, soliddiamond], numpoints=9, adaptive=false):

C:=plots:-textplot([[0.05,-1,F=0.1],[0.05,1.5,F=0.2], [0.05,3,F=0.5]], font=[times, 14]):

plots:-display(A, B, C, size=[500,500]);
```

## OK...

@vv  Formally, you are right. Because the velocity of motion in this model is equal to the gradient, we obtain a solution that is formally equal to  (0,0)  only for  t=infinity. But for  t=10  we get  (x,y) = (10*exp(-40),10*exp(-20)) = (4.248354255*10^(-17), 2.061153622*10^(-8)) , so that you are already at the top.

It is easy to write equations where the speed is constant, but this complicates the model.

## ?...

I did not understand the meaning of this

## Re...

@digerdiga  Do you want this output?

ex:=exp(-x^2/2*epsilon)*%sqrt((sqrt(x)/2/h)^2);