Carl Love

## 26715 Reputation

11 years, 249 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## truncate the sum...

Maple is unable come up with a numeric approximation of the infinite sum. (Maple's ability with numeric infinite summation is much more limited than its ability with numeric integration.)

Your sum's nth term is of the form O(1)*exp(-A*n^2*t), where A is a positive constant. (Don't worry if you don't know what O(1) means.) The n^2 in the exponent makes the terms get very small very fast (assuming t > 0). So, those terms with, say, n > 20 make no appreciable difference in the sum.  So, just change the sum's upper limit from infinity to 20. The 20 is somewhat arbitrarily chosen, but I couldn't find any difference at 15 digits between 20 and 21. That's far more accuracy than you need for a plot. If you need a concrete proof that truncating the sum at a certain n makes no difference at a certain number of digits, let me know.

I don't think a logplot is appropriate for this function because the range of function values is quite small. Also, if you don't provide a range for the independent variable, t, the plot commands will use -10..10. So, you want something like

plot(eq, t= 0..10^4);

Remember that the truncation is only valid for t > 0.

I'm not completely confident with my solution given below. So I'd appreciate some commentary from some of the other experts.

 > restart:
 > eqs:= {      R11+R12+fx = 1, R21+R22+fy = 1, R31+R32+fz = 1,      R11+R21+R31 = 1, R12+R22+R32 = 1,      R11*R22*fz + R12*fy*R31 + fx*R21*R32 -           R31*R22*fx - R32*fy*R11 - fz*R21*R12      = 1 }:
 > vars:= {R||(1..3)||(1..2)};

 > solve(eqs, vars);

No solution returned. The command eliminate is a lot like solve, but may provide some additional information.

 > Sol:= eliminate(eqs, vars);

Note the derived residual equation: fx + fy + fz = 1 (at the end). That's why solve didn't return anything. Also note that one of the Rs is left free. Let's put R11 into the basis (If one of the Rs has to be free, it might as well be one of the ones that we need to minimize.) and take fz out of the basis (Taking out one of the fs will cause eliminate to incorporate the residual equation into the solution.).

 > Sol:= eliminate(eqs, vars union {fz} minus {R11});

The expression to minimize:

 > Tr:= eval(R11+R22+fz, Sol[1]);

The expression is linear in R11; hence its minimum must occur at a boundary value for R11.

## print and colon...

Suppress the output of the loop by ending with a colon (:) instead of a semicolon (;). Unfortunately, this can only be applied to the whole loop; surprisingly, it makes no difference what punctuation you end the solve with. The colon after the loop will suppress all of its output, so now you force the output that you want by using print.

for alpha from .4 by .1 to 5 do
S := solve(x^4-alpha, x);
print(S[2]) ;
end do:

## gridrefine, crossingrefine...

The smartplot is using a finer mesh or some other refinement. If you add the option gridrefine= 2 (or higher) or crossingrefine= 2 (or higher) to your implicitplot, you will get a much better plot. While these may seem like expensive options, I point out that the plot below takes 47 millisecs for me, and the plot structure only stores 104 points.

`plots:-implicitplot(     (x^2+y^2)^2 = x^3+3/10*x*y^2, x= 0..1, y= -1..1,     gridrefine= 3, crossingrefine= 3,      resolution= 2^11, adaptranges,     scaling= constrained);`

## print(plot1)...

After the command you have above, put the command

print(plot1);

I think that that's all you need to do. Your command looks correct. But, not having the worksheet, I can't be sure of that.

## Better command is isolve...

The command ?isolve is better than msolve for this situation. But I don't think that there is any solution. Do long division. First factor 3 from the denominator:

3*p = (k^2-87)/(k+39) = k - 39 + 2*3*239/(k+39)

That leaves only a small number of possibilities for k, and I checked them all.

## With an Array-initializing procedure...

One of many ways to do it is with an Array-initializing procedure:

`curve:= dsolve(     {diff(y(x),x) = 1/(x-4.98), y(0)=1},     numeric, output= listprocedure):(Note that I used listprocedure.)pts:= < map(     P-> Array(          0..50,           proc(m) try P(m/10) catch: undefined end try end proc,          datatype= float[8]     ),     map(rhs, curve))[]>^%T:The resulting Array is plotable as is, because plot doesn't care if a few of your points are undefined; it just ignores them.plot(pts);`

## Next 55 terms...

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

## Save the value of omega on each iteratio...

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 ifend do; `

## Spherical coordinates...

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).

## 0...

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.

## Variance = (Standard deviation)^2...

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.

## Differ by a constant...

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.

## parameterized colours with colour= HUE(....

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.`

## both sides?...

I don't know what you mean by "both sides". Do you mean to solve it for each variable, T0 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});

>