Carl Love

Carl Love

26832 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The above Answers assume that the integral is real as if that were obvious. Perhaps it is obvious to some, but it wasn't to me. So here's a proof that it is real.

 

J:= exp(I*epsilon*ln((2+t)/(2-t)))*exp(I*n*t)/sqrt(4-t^2)/Pi;

exp(I*epsilon*ln((2+t)/(2-t)))*exp(I*n*t)/((-t^2+4)^(1/2)*Pi)

plot(eval(Im(J), [epsilon= -0.4, n= 1]), t= -2..2);

The imaginary part appears to be an odd function of t. That would make the integral over a symmetric real interval real. Let's prove it's an odd function.

evalc(Im(J)) assuming t::RealRange(-2,2);

(sin(epsilon*ln((2+t)/(2-t)))*cos(n*t)+cos(epsilon*ln((2+t)/(2-t)))*sin(n*t))/((-t^2+4)^(1/2)*Pi)

combine(%);

sin(epsilon*ln(-(2+t)/(t-2))+n*t)/((-t^2+4)^(1/2)*Pi)

ImJ:= unapply(%, t):

simplify(ImJ(t)+ImJ(-t)) assuming epsilon::real, n::real, t::RealRange(-2,2);

0

Therefore the imaginary part of the integrand is an odd function of t (for real n and epsilon).

 

Download odd_imag_part.mw

Why use such a crude trial-and-error method to find the B for each A? Here's a way using Optimization:-Minimize:

restart:
for k to 300 do
     A:= .01*k;
     Bmin:= Optimization:-Minimize(
          B, {exp(B)-A-1 >= 0, exp(B)-A/2-1.5 >= 0},
          feasibilitytolerance= 1e-6, optimalitytolerance= 1e-6,
          initialpoint= {B= .1}
     )[1];
     Sol[k]:= [A,Bmin]
end do:
plot(convert(Sol,list), labels= ['A', B[min]]);

Or we can get a fully symbolic solution thus:
L1:= solve(exp(B)-A-1, B);
                           ln(A + 1)
L2:= solve(exp(B)-A/2-3/2, B);
                            /1     3\
                          ln|- A + -|
                            \2     2/
convert(max(L1,L2), piecewise) assuming A>0;
                 /          /1     3\                  \
        piecewise|A <= 1, ln|- A + -|, 1 < A, ln(A + 1)|
                 \          \2     2/                  /


We can give this problem a more sophisticated combinatorial treatment. This might be appropriate for a case where the number of permutations is so great that iterating through them is undesirable. We note that any permutation of the four fractions in any solution is also a solution, essentially the same. Thus we can reduce the search by a factor of 4! = 24. The following returns the 8 unique solutions ((Kitonum's 192) / 4!) in 0.343 secs.

restart:
eq:= add(a[k]/a[k+4], k= 1..4) = a[9]:
S:= {$1..9}:
Sols:= table():

for s in S do
     S1:= S minus {s};
     for C in combinat:-choose(S1,4) do
          C1:= S1 minus C;
          for P in combinat:-permute([C1[]]) do
               A:= [C[],P[],s];
               if evalb(eval(eq, a= A)) then
                    Sols[A]:= [][]
               end if
          end do
     end do
end do:
 
for Sol in indices(Sols, nolist) do
   print(eval(eq, a= ``~(Sol)))
end do;

                  (2)   (3)   (5)   (7)      
                  --- + --- + --- + --- = (9)
                  (8)   (6)   (4)   (1) 

     
                  (1)   (4)   (6)   (7)      
                  --- + --- + --- + --- = (5)
                  (3)   (8)   (9)   (2)  

    
                  (3)   (5)   (6)   (7)      
                  --- + --- + --- + --- = (9)
                  (2)   (1)   (8)   (4)   

   
                  (3)   (4)   (5)   (7)      
                  --- + --- + --- + --- = (8)
                  (9)   (1)   (2)   (6)   

   
                  (2)   (4)   (5)   (7)      
                  --- + --- + --- + --- = (9)
                  (3)   (8)   (6)   (1) 

     
                  (2)   (5)   (6)   (9)      
                  --- + --- + --- + --- = (7)
                  (1)   (4)   (8)   (3)  

    
                  (2)   (5)   (6)   (7)      
                  --- + --- + --- + --- = (9)
                  (8)   (1)   (3)   (4)    

 
                  (1)   (2)   (7)   (9)      
                  --- + --- + --- + --- = (5)
                  (6)   (8)   (3)   (4)      

 

First, get rid of with(linalg). Then replace these old linalg commands

with the single command

Solution:= LinearSolve(A, XY):

Then you can plot the exact and the approximate solution together with

F:= x-> (1-x)*cos(x):
display([
     pointplot([seq]([h*(k-2), Solution[k]], k= 2..m+1), color= blue),
     plot(F(x), x= 0..1)
], symbol= diamond
);

You can plot the errors with

pointplot(
     [seq]([h*(k-2), F(h*(k-2)) - Solution[k]], k= 2..m+1),
     symbol= cross, title= "Errors"
);

You need to use lists rather than sets because you have no control over the order that things will appear in a set. Also, your indices themselves need to be lists because "naked" sequences cannot be distinct members of lists (or sets). Once you have lists, let's say A is your list of indices and B is the set of values. Then just do

A =~ B

It is possible for a naked sequence to be either or both of the sides of an equation. You can achieve this by

zip((x,y)-> op(x)=y, A, B)

zip is the two-set (or two-list) analog of map that you were looking for.

Example:

A:= [[1,2], [3,4]]:
B:= [5,6]:
zip((x,y)-> op(x)=y, A,B);
                    [(1, 2) = 5, (3, 4) = 6]

There is no effective way to make such an assumption. Even if you were include every possible inequation x[i] <> x[j], you would likely just overload the ability of a solver. That's why Combinatorial Optimization and Constraint Satisfaction Problems are their own areas of study, distinct from, for example, Linear Programming.

But you can solve such an equation by cycling through all permutations if the number is not too great. In the case at hand, we have 9! ~ 300,000---reasonable. The following executes in 0.235 secs for me.

restart:
eq:= a[9] = add(a[2*k-1]/a[2*k], k= 1..4):
A:= combinat:-firstperm(9):
while not evalb(eval(eq, a= A)) do  A:= combinat:-nextperm(A)  od:
eval(eq, a= ``~(A));

I tried in both Maple 16 and 17, and I cannot duplicate your results; I get relatively smooth plots for the first and second derivative. According to ?dsolve,numeric,bvp , the maximum value of maxmesh is 2^13 = 8192. I got results at 2^10, with no difference at 2^13. I decreased abserr to 1e-7 (default is 1e-6) and increased Digits to 20, and I got no difference. Here are my plots. Note that the range is 0..1.

for k from 0 to 2 do
     odeplot(sol, [x, diff(Urn(x), [x$k])], 0..1, numpoints= 1000)
end do;

You wrote:

Also, the list with the output does not show the 2nd derivative of Uthetan, the 4th derivative of Urn and the 2nd derivative of Uxn. They are in the equations... is that a problem?

It is not a problem. It never includes the highest-order derivative of any of the functions. If you want to plot those derivatives, there are several workarounds.

I tried in both Maple 16 & 17, obtaining the same results. The "hang" occurs in simplify. I don't know whether it is an infinite loop, or just an expression that is too complicated for simplify. Since the expression fits on one line on my screen, the latter seems unlikely. Running with infolevel[simplify]:= 5, the only information is a seemingly endless repetition of

  • simplify/do: applying  commonpow  function to expression
  • simplify/do: applying   power  function to expression

The hang will occur for any of z[17], z[18], z[19], z[20]. It does not matter whether the expand or factor or both are applied before the simplify. Only the first two solutions need to be substituted in for the hang to occur. So here is one of the seemingly simple expressions for which it hangs:

subs(solution1, solution2, z[20]);

 

lprint(%);


0 = -(-b4*f8+b8*f4)*conjugate(f7)*d5*(conjugate(b4)*conjugate(d8)-conjugate(b8)*conjugate(d4))/((a4*b8-a8*b4)*(conjugate(a4*b8)-conjugate(a8*b4)))-(-a4*f8+a8*f4)*conjugate(f7)*d5*(conjugate(a4)*conjugate(d8)-conjugate(a8)*conjugate(d4))/((a4*b8-a8*b4)*(conjugate(a4*b8)-conjugate(a8*b4)))+conjugate(e7)*e5

 

 

How about diff(z(t), t) = m with initial condition z(0) = n? Then the system becomes

diff(x(t), t) = -x(t) + 2*z(t)*y(t) - 2*k*z(t)^2,

diff(y(t), t) = -y(t) + k*m - z(t)*(x(t)/2 - 1),

diff(z(t), t) = m;

Here it is, by guiding Maple at each step (or "holding its hand").



 

restart:

J:= n-> int(sec(x)^n, x);

proc (n) options operator, arrow; int(sec(x)^n, x) end proc

IntegrationTools:-Parts(J(n), sec(x)^(n-2));

sin(x)*sec(x)^(n-2)/cos(x)-(int(sin(x)*sec(x)^(n-2)*(n-2)*tan(x)/cos(x), x))

algsubs(sin(x)/cos(x)= tan(x), %);

tan(x)*sec(x)^(n-2)-(int(sec(x)^(n-2)*(n-2)*tan(x)^2, x))

algsubs(tan(x)^2= sec(x)^2 - 1, %);

tan(x)*sec(x)^(n-2)-(int(sec(x)^(n-2)*(n-2)*(sec(x)^2-1), x))

IntegrationTools:-Expand(%);

tan(x)*sec(x)^(n-2)-n*(int(sec(x)^n, x))+n*(int(sec(x)^n/sec(x)^2, x))+2*(int(sec(x)^n, x))-2*(int(sec(x)^n/sec(x)^2, x))

combine(%, power);

tan(x)*sec(x)^(n-2)-n*(int(sec(x)^n, x))+n*(int(sec(x)^(n-2), x))+2*(int(sec(x)^n, x))-2*(int(sec(x)^(n-2), x))

solve(J(n) = %, {J(n)})[];

int(sec(x)^n, x) = (tan(x)*sec(x)^(n-2)+n*(int(sec(x)^(n-2), x))-2*(int(sec(x)^(n-2), x)))/(n-1)

collect(%, J(n-2));

int(sec(x)^n, x) = (n-2)*(int(sec(x)^(n-2), x))/(n-1)+tan(x)*sec(x)^(n-2)/(n-1)

 

 



Download int_sec^n.mw

 

 

 

You need to create an Array (or Matrix) before you can assign to it. Otherwise, the object becomes a table.

(n,m):= (3,3);

U:= Array(0..n+2, 0..m):

U[.., 0]:= < y(a), seq(f[k], k= 0..n), y(b) >;

Note that it is U[.., 0] rather than U[0, ..]. The latter would create the first row rather than the first column.

First, you must create and store S as a Vector, list, or table---not a set. There are two problems with using a set:

  1. Duplicate entries are eliminated, and in your case there is a small chance that two of the numbers associated with an index pair are the same.
  2. A set is stored sorted, and you'll want to avoid the super-linear overhead of the sort.

Which one is the most efficient to use---Vector, list, or table---depends entirely on the code used to create S. So, I'd like to see that code. Note that if you use a table, then the sum can be accumulated as the entries are created. Create it with table(sparse) so that entries are initialized to 0.

The process of converting to a Matrix is independent of which of the three formats that you use. I will assume that the Matrix order N is known a priori; if instead the order needs to be determined from the scan of S, then my code will need to be modified.

 

restart:

Accumulate:= proc(L::{list, rtable, table})
     local e, T:= table(sparse);
     for e in L do  T[lhs(e)]:= T[lhs(e)] + rhs(e)  end do;
     T
end proc:

Example of its use:

n:= 1:  N:= 2^n:

R:= op@RandomTools:-Generate(list(posint(range= N), 2),

            makeproc

       )
    = RandomTools:-Generate(float(method= uniform), makeproc)

:

S:= ['R()'$ 2*N^2];

[(1, 1) = .661561962545309, (2, 2) = .244155091492448, (1, 1) = .632324605279485, (1, 2) = .881385194785638, (2, 2) = .586801202666963, (1, 1) = .137750180938911, (2, 2) = .250797100906032, (1, 1) = .891461629197522]

Matrix(N, Accumulate(S), storage= sparse, datatype= float[8]);

Matrix(2, 2, {(1, 1) = 2.32309837796122, (1, 2) = .881385194785638, (2, 1) = 0, (2, 2) = 1.08175339506544})

I will try to verify experimentally that this process's time complexity is O(|S|) and post a followup.

 

 

Download sparse.mw

Because of the issue of the target lists possibly appearing as sublists, I recommend that you use strings instead of lists as your data structure. Then you can use StringTools:-SubstituteAll for the operation.

Even if you need to use lists, it may be faster to perform the operation by converting to a string and then back to a list every time. That's because the StringTools commands are written in compiled code.

First, get rid of linalg. If you need linear algebra (you don't for this bit of code), use LinearAlgebra instead.

Only the terminal punctuation of the whole for-do loop, either colon or semicolon, makes any difference in the printed output. So we need to suppress all the output by ending the loop with a colon, and, within the loop, selectively printing what we want with a print statement.

for h from .5 by .1 to 1 do
     eq1 := x*(-b*x^2-x+1);
     eq2 := y*((a*x*x)/(b*y^2)-d-h*y);
     S := solve({eq1, eq2}, {x, y});
     SS := solve(subs(S[3], {omega^4+(h*y+x)*omega^2+h^2*x-y}), {omega});
     tau := simplify(subs(S[3], subs(SS[3], (b^2*h*y+a*x)/omega)));
     print('tau' = tau);
end do:


You wrote:

what if we change the original equation to a pde with the derivative on the left hand side to a derivative with repect to t?

Then the rest of your code in your original question works. The animation shows a very interesting wave pattern. Thus I think this is what the equation was origianlly intended to be. It doesn't make sense to me to do it in 3D though.

S:= pdsolve(pde, {bc, ic}, numeric):
S:-animate(t= 0.1, frames= 64);

First 355 356 357 358 359 360 361 Last Page 357 of 385