I don't have a formula. I am guessing based on a mental analysis of the pattern. No polynomials or computerized interpolation used. The next 55 terms are

59051, 59053, 59055, 59059, 59061, 59067, 59077, 59079, 59085,
59103, 59131, 59133, 59139, 59157, 59211, 59293, 59295, 59301,
59319, 59373, 59535, 59779, 59781, 59787, 59805, 59859, 60021,
60507, 61237, 61239, 61245, 61263, 61317, 61479, 61965, 63423,
65611, 65613, 65619, 65637, 65691, 65853, 66339, 67797, 72171,
78733, 78735, 78741, 78759, 78813, 78975, 79461, 80919, 85293,
98415

You just need to save the value of omega in the same way that you save KR. I put lines of ### around the code below that I changed. Let me know if this solves your problem.

Unrelated to your problem: There's no need to use evalm. The command serves no purpose in contemporary Maple.

In the future, please post your code without the ">" beginning each line. I have to remove those one by one to use your code.

Digits:=30;

 h0:=0.156;
 d:=0.32*h0;
 l:=2;
 h1:=h0-d;
 h2:=h0+d;
 h3:=0.6*h0;
 g:=9.8;
 d1:=1;
 Term:=8;
 Num:=100:
 n:=2:
 l1:=l/n;
 
 for N from 1 to Num do
 lambda:=2*n*Pi/l: 
 k0:=evalf(10^(-10)*Pi+2.5*(N-1)*Pi/(Num-1)):
 tau0:=evalf(k0*h0):
 omega:=evalf((g*k0*tanh(k0*h0))^(1/2)):
 ###############
 Omega[N]:= omega:
 ################
 E:=g/(omega)^2:
 k1:=abs(fsolve(omega^2=g*k*tanh(k*h1),k)):
 tau1:=evalf(k1*h1):
 k2:=abs(fsolve(omega^2=g*k*tanh(k*h2),k)):
 tau2:=evalf(k2*h2):
 k3:=abs(fsolve(omega^2=g*k*tanh(k*h3),k)):
tau3:=evalf(k3*h3):
 
 u0:=evalf(g*tanh(tau0)*(1+2*tau0/(sinh(2*tau0)))/(2*k0));
 u01:=evalf(g*tanh(tau3)*(1+2*tau3/(sinh(2*tau3)))/(2*k3));
 u1:=evalf(g*(1-(tanh(tau0))^2)*(sinh(2*tau0)-2*tau0*cosh(2*tau0))/(4*(2*tau0+sinh(2*tau0))));
 H:=evalf((1+2*tau0/sinh(2*tau0))/(-lambda*k0*d));
 delta00:=evalf((lambda*d*u1/u0+I*k0)*H);
 delta01:=evalf((lambda*d*u1/u0-I*k0)*H);
 delta11:=evalf((lambda*d*u1/u0+I*k0)*H*exp(I*k0*l));
 delta12:=evalf((lambda*d*u1/u0-I*k0)*H*exp(-I*k0*l));
 delta21:=evalf(exp(I*k0*(l+l1)));
 delta22:=evalf(u0*k0*exp(I*k0*(l+l1)));
 delta31:=evalf(exp(-I*k0*(l+l1)));
 delta32:=evalf(-u0*k0*exp(-I*k0*(l+l1)));
 delta41:=evalf(exp(I*k3*(l+l1)));
 delta42:=evalf(u01*k3*exp(I*k3*(l+l1)));
 delta51:=evalf(exp(-I*k3*(l+l1)));
 delta52:=evalf(-u01*k3*exp(-I*k3*(l+l1)));
 delta61:=evalf(exp(I*k3*(l+l1+d1)));
 delta62:=evalf(u01*k3*exp(I*k3*(l+l1+d1)));
 delta71:=evalf(exp(-I*k3*(l+l1+d1)));
 delta72:=evalf(-u01*k3*exp(-I*k3*(l+l1+d1)));
 delta81:=evalf(exp(I*k0*(l+l1+d1)));
delta82:=evalf(u0*k0*exp(I*k0*(l+l1+d1))); 
 
Hi:=(Matrix([[1,1],[delta00,delta01]]))^(-1);
H1:=(Matrix([[delta21,delta31],[delta22,delta32]]))^(-1);
H2:=Matrix([[delta41,delta51],[delta42,delta52]]);
H3:=(Matrix([[delta61,delta71],[delta62,delta72]]))^(-1);
H4:=Matrix([[delta81],[delta82]]);
 
########## 
MM:= H1.H2.H3.H4:
##########
 R:=evalf(MM[2,1]/MM[1,1]):
Kr:=abs(evalf(R)):
KR[N]:=abs(Kr);
 
 end do:
 
########### L:=seq(KR[N],N=1..Num):
 
############
for N from 1 to Num do
    if KR[N]>0.2 then 
        print(Omega[N])
        ###################
    end if
end do;
 

The integral is too difficult for Maple in its current form. Convert it to spherical coordinates.

I say this with trepidation, considering how badly I erred with the limits on your previous question: Are you sure that that set R is correctly specifed? The first condition, z <= sqrt(x^2+y^2+z^2), is satisfied by any three real numbers. The second condition is strange because there is no lower limit to z. I wonder if it was supposed to be abs(z).

The integral is easily seen to be 0 without the help of Maple.

2*y >= x^2 + y^2  -->  0 >= x^2 - 2*y + y^2 = (x-y)^2  -->  (x-y)^2 = 0  -->  x=y.

Thus D is meager, and the integral is 0.

Your code contains the expression

Variance(X)=sqrt(StandardDeviation(A)).

That should be

Variance(X) = StandardDeviation(A)^2.

Also, consider using ?KernelDensity rather than EmpiricalDistribution, which gives a discrete distribution. Then you can make a custom ?Distribution by uing the function returned by KernelDensity as the ?PDF.

The two answers differ by a constant, so they are both valid forms of the integral. The constant is ln(-1) = Pi*I. You can take the derivative of each and see that in both cases you get the integrand.

It is not a silly question. The easiest way to get parameterized colours is to use colour= HUE(x) where 1 >= x >= 0. The HUE scale is the standard colour spectrum, from red=0 to violet=1.

Here's an example. I've generated example data such as you describe, but using random data since I don't have your specific vectors.

V:= unapply(Vector(16, ()-> randpoly(x, degree= 1)), x):
V||(1..4):= 'V(k)' $ 'k'= 1..4:
So now I have 4 vectors, V1, ..., V4 of 16 numeric elements each.
plot([seq]([seq]([k, (V||k)[j]], k= 1..4), j= 1..16), colour= [seq](HUE(j/16), j= 1..16));
Note that the argument to HUE is matched to the index into the Vectors.

I don't know what you mean by "both sides". Do you mean to solve it for each variable, T₀ and t? That can be done like this (I've changed the decimals to exact fractions, but it is not necessary to do that):

> 



eq:= (T[0]+exp(-sqrt(t/1000)))/T[0] = 95/100;

















> 



solve(eq, {T[0]});

















> 



solve(eq, {t});

















> 



 


















Download solve2.mw

Working from the Wikipedia pages "Gauss-Hermite quadrature" and "Hermite polynomials", I whipped up this procedure for you. It takes a positive integer argument n and returns two lists of n elements, the first being the nodes and the second being the weights.

GaussHermite:= proc(n::posint)
local 
     x,
     Hm:= unapply(numapprox:-hornerform(orthopoly[H](n-1, x), x), x), 
     R:= [fsolve](orthopoly[H](n,x), x),
     C:= evalf(2^(n-1)*n!*sqrt(Pi)/n^2)
;
     R, map(x-> C/Hm(x)^2, R)
end proc:
     
GaussHermite(6);
 [-2.35060497367449, -1.33584907401370, -0.436077411927616, 
   0.436077411927616, 1.33584907401370, 2.35060497367449], 

   [0.00453000990550894, 0.157067320322850, 0.724629595224395, 
   0.724629595224395, 0.157067320322850, 0.00453000990550894]

Your question is not boring, and you've done nothing wrong. You've stumbled upon an obscure bug in Maple's arbitrary-order derivatives.

Are you really trying to "solve" the formula in some way? Or do you just want to evaluate it for specific integer values of r? If the latter:

Maple tries to take the derivative of arbitrary order that occurs in your sum. The formula that it returns is apparently wrong (very wrong, in fact). We can work around that by using add instead of sum.

Mrs:= proc(s, r::posint, x)
local 
     i, X,
     S:= s+1,
     P:= i-> i+s,
     T:= i-> 2*i+s,
     U:= i-> 2*i+s+1,
     R:= t-> (2*t)!/t!/4^t,
     A:= R@P, B:= R@T, C:= R@U
;
     eval(
          add(
               (-1)^i*(A/B*C)(i)/T(i)!/i!/4^i*
               diff((1-X^2)^T(i), [X $ 2*i]),
               i= 0..r-1 
          ),
          X= x
     )
end proc:

If you are trying to "solve" the formula in some way other than simple evaluation, let me know. My procedure does allow you to use symbolic or non-integer s, and some significant and interesting simplifications are possible (try s = 1/2).

You wrote:

While trying your solution, I noticed that my maple gives me a higher order of the complex part of the solution (E-6).

That's because I had Digits set to 15 and you had Digits set to 10. You can include the command Digits:= 15 right after your restart.

You wrote:

Furthermore I cannot plot my solution due to the following error....

Your plot command has unbalanced quotation marks, and possibly some other bad character(s). I'm surprised that it got past the parser. I recommend that you retype (not cut and paste) the whole command in a new execution group. All of the quotes need to be double quotes ("), not single quotes ('). I also recommend that you use input mode Maple Input (like the red characters that I used) rather than 2D Input, because with Maple Input both you and I can see exactly what you typed. With 2D Input it is impossible to visually distinguish a double quote from two single quotes next to each other.

You wrote:

Furthermore I get an error trying to find the slope of the catenary at x=0.

Use eval(Helling, x= 0).

You wrote:

I was wondering why maple supplies us with two solutions for the differential equation.

There is only one solution to the differential equation proper; it is only after the boundary conditions are put in that the branching occurs. You can see that the differential equation has one rather simple solution by calling dsolve without the boundary conditions:

dsolve(DIV);
                       cosh(_C1 qT + qT x)      
                y(x) = ------------------- + _C2
                               qT               

You wrote:

Is it because of the square root in the system?

No. When you apply the boundary conditons to the above solution to solve for the constants of integration _C1 and _C2, you get an equation of the form cosh(x) = B. Such equations generally have two solutions (for B > 0) because horizontal lines cross cosh curves twice.

You wrote:

Furthermore, what is the origin of the complex part of the solution?

After the boundary conditions are applied, the solution contains an expression of the form cosh(A+ln(RootOf(Q))) where Q is a quadratic with two real roots, one positive and one negative. (The nature of the roots depends on your parameter values and boundary conditions. For the case we are discussing, one is positive and one is negative.) The ln of a negative number is complex; specifically, for x < 0, ln(x) = ln(|x|) + Pi*I.

It is important to note that this is not a complex solution; we are not ignoring any part of the true solution. The Pi*I turns real after the cosh is applied to it. It is only because of the inexact floating-pointing arithmetic that we see it at all. That's why I call them spurious imaginary parts. Below, I redid the worksheet in such a way that exact arithmetic is used up to the point that the complex parts disappear.

Download SpuriousImag2.mw

@emma hassan

Your last command has a product of two vectors, each with five elements (assuming m = 4), which you are trying to assign to a vector. But the product of two vectors ordinarily means the dot product, which yields a scalar (a single number, not a vector). Did you mean for that to be the elementwise product of the two vectors? The command for that would be

U[j+1, ..]:= U[j, ..] *~ A;

(assuming that you're using Maple 13 or later. If you have an earlier Maple, let me know.)

Also, the loop should go for j from 0 to n-1, not for j from 0 to n.

You asked that, and it was answered, within the past week. If ex is your expression, then

series(ex, c= infinity, 3);

Did that answer not work for you?

Thanks for using a shorter title.

You cannot use an index of 0 in a Maple Matrix. Matrix indices always start at 1. But you can use 0 indices in the closely related structure called Array.

h:= 1/4:
m,n:= 4,4:
U:= Array(0..m, 0..n, datatype= float[8]):
U[0, ..]:= < seq(evalhf(cos(Pi*i*h)), i= 0..m) >;

Note that I used evalhf instead of evalf. The usual evalf would work also, but I think that evalhf is a little better in this situation.