mmcdara

3713 Reputation

17 Badges

5 years, 347 days

MaplePrimes Activity


These are answers submitted by mmcdara

@VK057 

Here are two methods, probably amid several others:

restart:

f := 2.0*cos(1.5*x + 4.0) + 2.2:

g := 2.0*cos(1.6*x + 3.85) + 2.0:

 

# Method 1

r1 := fsolve(f=g, x=-1..0);
r2 := fsolve(f=g, x= 0..4);

 

-.7967442905

 

3.834531622

(1)

# Method 2

eq := unapply(f-g, x);
r1 := RootFinding:-NextZero(x -> eq(x), -1);
r2 := RootFinding:-NextZero(x -> eq(x), r1);

proc (x) options operator, arrow; 2.0*cos(1.5*x+4.0)+.2-2.0*cos(1.6*x+3.85) end proc

 

-.7967442905

 

3.834531622

(2)

Download TwoMethods.mw

More generaly, if you want to select all the roots of f-g within some interval [x1, x2] (here [-8, 8]) you can do this

x1 := -8:
x2 :=  8:
r  := NULL:
while x1 < x2 do
  x1 := RootFinding:-NextZero(x -> eq(x), x1);
  r  := r, x1
end do:
r := select(`<`, [r], x2)
   [-6.502888690, -4.665889115, -2.362208506, -0.7967442905, 

     3.834531622, 5.406998312, 7.708001175]

Consider this loop

for i from 0 to N do:  
dsolve({equ[1][i],cond[1][i]},X[i](t));  
X[i](t):=rhs(%):  
end do;

and look to the ode you solve:

i:=1:
{equ[1][i],cond[1][i]}
      // d         \                                       
     { |--- X[1](t)| + 2 X[0](t) Y[0](t) = 0, X[1](0) = 0, 
      \\ dt        /                                       

                                                  \ 
       X[1](5) = 0, D(X[1])(0) = 0, D(X[1])(5) = 0 }
                                                  / 

This is a first order ode with 4 boundary conditions while only one should be given.
More of this the boundary condition cannot be a condition on the derivative of X[1](t).
As a result dsolve returns "nothing" and the command

rhs(%)

generates an error.

It is up to you to know what is the correct boundary condition 

X[1](0) = 0;
#or
X[1](5) = 0;

 

First point: 
A simple way to generate your system of equations is:
 

N   := 4:   # or 9 in your last reply
P   := 2*N:
sys := [ add(W__||n, n=1..N)=4, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
print~(sys):  # check


Second point;
The coefficients 2/3, 2/5, ..., 2/(2*p+1) are (curiously ?) the values of the integral of the pth Legendre polynomial squared over -1 to 1:

seq([p, int(simplify(LegendreP(p, x),'LegendreP')^2, x=-1..1)], p=1..6)
        [   2]  [   2]  [   2]  [   2]  [   2 ]  [   2 ]
        [1, -], [2, -], [3, -], [4, -], [5, --], [6, --]
        [   3]  [   5]  [   7]  [   9]  [   11]  [   13]

Your set of equations has an an extremely strong resemblance (apart for the first relation W__1+...+W__N=4 instead of the "normal" W__1+...+W__N=2 relation) with the set of relations that enable computing the weights (W) and the locations (zeta) of the Gauss-Legendre points.
Is this a coincidental resemblance ?

For instance this code returns the weights and locations of the Gauss-Legendre points for a N points Gauss quadrature.

N   := any_strictly_positive_integer:
P   := 2*N:
sys := [ add(W__||n, n=1..N)=2, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
print~(sys): #check
solve(sys);


These same points can be computed in a far more efficient way for large values of N
(see for instance Gauss–Legendre_quadrature or 9403469.pdf)

Thus, if your  "first" equation W__1+...+W__N=4 is indeed the correct one, one can imagine that the methods mentioned in this link could be adjusted to solve your specific system.


Third point;
If the "idea" suggested at the end of the previous point is a dead end, one can simplify your problem by assuming that there exist symmetries of "Gauss-Legendre points".
For instance, without loss of generality, and if N=4, on can assume that 

W__1 = W__4, W__2 = W__3, zeta__1 = -zeta__4, zeta__2 = -zeta__3

Next one can infer that W > 0 and that abs(zeta) <=1 (for the rhs of relations of odd ranks are decreasing as the powers of zeta are increasing).

For instance:

N   := 6:
P   := 2*N:
sys := [ add(W__||n, n=1..N)=N, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:

t0  := time():
sol := fsolve(sys): # no solution found
time()-t0;
                             16.204

#----------------------------------------------
hypothesis := [seq(W__||(N-i+1)=W__||i, i=1..N/2), seq(zeta__||(N-i+1)=-zeta__||i, i=1..N/2)]:
redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):

t0  := time():
redsol := fsolve(redsys): # no solution found
time()-t0;

                             2.753
#----------------------------------------------
t0  := time():
redsol2 := fsolve(redsys, {seq(W__||i=0..4, i=1..N/2), seq(zeta__||i=-1..1, i=1..N/2)});
time()-t0;
{W__1 = 2.432885814, W__2 = 0.1793486330, W__3 = 0.3877655527, 

  zeta__1 = -0.08695232853, zeta__2 = -0.9294284522, 

  zeta__3 = -0.6423774108}
                             0.025


NO MIRACLE HERE: for N>6 no solution is found and some work remains to be done
For instance :

  1. Transform the problem into a minimization one.
  2. When N is odd, forcing the "central" zeta to 0  (zeta__i where i=ceil(N/2))
  3. Taking into account that for the values of zeta for N-1 and N are alternate in the sense that
    • -1 < zeta__1 < zeta[N-1, 1]
    • zeta[N-1, , k] < zeta__k < zeta[N-1, k+1] for k=2..N-1
    • zeta[N-1, N-1] < zeta__N < 1

This could provide better search ranges for the zeta

Here is procedure aimed to compute W and zeta for relatively large values of N which implements these ideas.
It returns:

  • the solution as a list of equalities of the form W__i =...  and  zeta__i = ...
  • the list of absolute errors when this solution is plugged into the system to solve.

The second argument in the Gale procedure is a switch to select the "natural" -1..1 range for zeta, or the "potentially optimal" ranges defined in point 3 above.
You will see in the attached file that the natural range option gives var better results.

GaLe.mw

restart

GaLe := proc(M, zoom)
  local N, P, sys, hypothesis, redsys, constr, sol, nm, xi, here, J, abserr:

  Digits := 20:
  N          := min(M, 6):
  P          := 2*N:
  sys        := [ add(W__||n, n=1..N)=4, seq(add(W__||n*(zeta__||n)^(p-1), n=1..N)=2/p*modp(p, 2), p=2..P) ]:
  hypothesis := [seq(W__||(N-i+1)=W__||i, i=1..N/2), seq(zeta__||(N-i+1)=-zeta__||i, i=1..N/2)];
  redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):
  if N::odd then
    redsys := eval(redsys, zeta__||(ceil(N/2))=0):
  end if:
  constr     := {seq(W__||i=0..4, i=1..ceil(N/2)), seq(zeta__||i=-1..1, i=1..N/2)}:
  sol        := fsolve(evalf(redsys), constr);

  if M > 6 then
    for nm from N+1 to M do
      xi   := sort(eval([seq(zeta__||i, i=1..nm/2)], sol)):
      here := [-1..xi[1], seq(xi[k]..xi[k+1], k=1..nops(xi)-1)];
      if nm::even then here := [here[], op(2, here[-1])..0]: end if:

      P          := 2*nm:
      sys        := [ add(W__||n, n=1..nm)=nm, seq(add(W__||n*(zeta__||n)^(p-1), n=1..nm)=2/p*modp(p, 2), p=2..P) ]:
      hypothesis := [seq(W__||(nm-i+1)=W__||i, i=1..nm/2), seq(zeta__||(nm-i+1)=-zeta__||i, i=1..nm/2)];
      redsys     := eval(sys[[seq(1..P, 2)]], hypothesis):

      if nm::odd then
        redsys := eval(redsys, zeta__||(ceil(nm/2))=0):
        hypothesis := [hypothesis[], zeta__||(ceil(nm/2))=0];
      end if:

      J      := evalf(add(((lhs-rhs)~(redsys))^~2)):
      if zoom then
        constr := {
                    seq(W__||i >= 0, i=1..ceil(nm/2)),
                    seq(W__||i <= 4, i=1..ceil(nm/2)),
                    seq(zeta__||i >= op(1, here[i]), i=1..floor(nm/2)),
                    seq(zeta__||i <= op(2, here[i]), i=1..floor(nm/2))
                  }:
      else
        constr := {
                    seq(W__||i >= 0, i=1..ceil(nm/2)),
                    seq(W__||i <= 4, i=1..ceil(nm/2)),
                    seq(zeta__||i >= -1, i=1..floor(nm/2)),
                    seq(zeta__||i <= +1, i=1..floor(nm/2))
                  }:
      end if:
      sol := Optimization:-Minimize(J, constr, iterationlimit=1000)[2];
      if nm::odd then sol := [sol[], zeta__||(ceil(nm/2))=0] end if:
    end do:
    sol := [sol[], eval(hypothesis, sol)[]];
  end if:  
  if M <= 6 then
    sol := [sol[], eval(hypothesis, sol)[], `if`(M::odd, zeta__||(ceil(M/2))=0, NULL)];
  end if:
  abserr := abs~(evalf(eval((lhs-rhs)~(sys), sol))):
  return sol, abserr:
end proc:

# 2nd argument of procedure Gale:
#    false means zeta is searched between -1 and +1
#    true  means zeta is searched in the interval defined by the values obtained for N-1

N := 9:
sol_F := CodeTools:-Usage( GaLe(N, false) ):

print():

sol_T := CodeTools:-Usage( GaLe(N, true) ):

memory used=348.60MiB, alloc change=36.00MiB, cpu time=2.03s, real time=4.70s, gc time=150.48ms

 

 

memory used=140.04MiB, alloc change=0 bytes, cpu time=811.00ms, real time=814.00ms, gc time=72.94ms

 

dataplot(
  <[sol_F[2], sol_T[2]]>,
  color=[blue, red],
  legend=["False", "True"],
  labels=["eq number", "abserr"],
  labeldirections=[default, vertical]
);

 

ix_F    := sort( eval([seq(zeta__||i, i=1..N)], sol_F[1]), output=permutation):
zeta_FS := [ seq(eval(zeta__||(ix_F[i-N]), sol_F[1]), i=N+1..2*N) ]:
W_FS    := sort(  eval([seq(W__||i, i=1..N)], sol_F[1]) ):
W_FS    := [ W_FS[[seq(1..N, 2)]][], W_FS[[seq(N..1, -2)]][] ]:

#ix_T := sort( eval([seq(zeta__||i, i=1..N)], sol_F[1]), output=permutation):

sol_TS := sort(sol_T[1]):

Matrix(
  2*N, 3,
  (i,j) -> `if`(
             j=1,
             `if`(i <= N, W__||i, zeta__||(i-N)),
             `if`(
               j=2,
               `if`(i <= N, W_FS[i], zeta_FS[N-i]),
               `if`(i <= N, eval(W__||i, sol_TS), eval(zeta__||(i-N), sol_TS))
             )
           )
):
evalf[6](%):
< `<|>`(["", "False", "True"]), % >:

# DocumentTools:-Tabulate(%, width=30, alignment=left):

plots:-display(
  seq(
    plot([[zeta_FS[i], -0.1], [zeta_FS[i], 0.1]], thickness=7, color=blue, transparency=0.7, legend="False"),
    i=1..N
  ),
  seq(
    plot([[eval(zeta__||i, sol_T[1]), -0.1], [eval(zeta__||i, sol_T[1]), 0.1]], thickness=3, color=red, legend="True"),
    i=1..N
  ),
  view=[-1..1, -0.1..0.1],
  axis[2]=[tickmarks=[]],
  scaling=constrained
)

 

# Legendre polynomials compared to your "pesudo-Legendre" polynomials

Matrix(
  3, 3,
  (m,n) ->
    plots:-dualaxisplot(
      plot(simplify(LegendreP(3*(m-1)+n, x),'LegendreP'), x=-1..1, color=blue, legend="Legendre"),
      plot(
        expand(eval(mul((x-zeta__||i), i=1..3*(m-1)+n), GaLe(3*(m-1)+n, true)[1])),
        x=-1..1, color=red, legend="pseudo Legendre"
      ),
      title=cat("N = ", 3*(m-1)+n)
    )
):

# DocumentTools:-Tabulate(%, width=60):

 

 

 

 

Export your plot, either by clicking on it or programatically (see plotsetup), in jpeg, png or pdf, and insert thisfile into your docx document

You want to solve a system of 4 linear equations depending on 7 parameters (a, b, c, t, x, y, z) by condidering it as a system of 4 linear equations in 3 unknowns (a, b, c), parameterized by (t, x, y, z).
A priori it's not surprising that you got no solutions.

Try any of these instead:
 

solve({seq(G[i] = 0, i = 1 .. 4)});
{a = 2 y + 2 z, b = -2 x + 2 y - c, c = c, t = x - y - z, x = x, 

  y = y, z = z}
for u in [t, x, y, z] do
solve({seq(G[i] = 0, i = 1 .. 4)}, {a, b, c, u});
end do;
   {a = 2 y + 2 z, b = -2 x + 2 y - c, c = c, t = x - y - z}
   {a = 2 y + 2 z, b = -c - 2 z - 2 t, c = c, x = y + z + t}
   {a = -2 t + 2 x, b = -c - 2 z - 2 t, c = c, y = x - z - t}
   {a = -2 t + 2 x, b = -2 x + 2 y - c, c = c, z = x - y - t}


Other solution, it it makes sense to you, solve your system in the "least squares sense" :

vars := [a, b, c]:
M,V := GenerateMatrix(sys, vars):
<vars> = (M^+.M).(M^+.V);

Just do this
 

# write the general solutions
Un := (x, t, n) -> T(t, n)*X(x, n);

# Derivative of Un w.r.t t
DtUn := eval(diff(Un(x, t, n), t), t=0):

# Compute the Fourier coefficients
int(DtUn*sin(Pi*n*x/l), x=0..l):
int(psi*sin(Pi*n*x/l), x=0..l) assuming l::real:
isolate(%%=%, C1[n]);
    
                                  2 /                 /1     \\
          piecewise(l < 0, 0, 1) l  |sin(Pi n) - 2 sin|- Pi n||
                                    \                 \2     //
C1[n] = - -----------------------------------------------------
                2  2 /  1                         1       \    
              Pi  n  |- - a cos(Pi n) sin(Pi n) + - a Pi n|    
                     \  2                         2       /    

 

If you realy want to use the notation e^(...) instead of the "normal" one exp(...) you can do this

restart:
local e := exp(1):
uu := (x, y, t) -> e^(.2*t)*(e^(-x)+e^(-y)):
uu(.1, .1, .1);
                               2      
                          ------------
                                  0.08
                          (exp(1))    
evalf(%)
                          1.846232693



Or, if you are only interested in floating point values 

restart:
local e := exp(1.):
uu := (x, y, t) -> e^(.2*t)*(e^(-x)+e^(-y));
uu(.1, .1, .1)
                          1.846232693

But listening acer is a far better idea

EDITED  10/27 9:30 pm GMT

For a reason I don't understand it seems that you must transform your sine expression before using solve.

restart:
f := sin(2*t)-sin(3*t)+sin(4*t);     
g := expand(f, trig):
sol_g := solve([g, `and`(t >= 0, t <= 2*Pi)], allsolutions, explicit=true)

 /    1   \    /    5   \    /    2   \    /    4   \            
{ t = - Pi }, { t = - Pi }, { t = - Pi }, { t = - Pi }, {t = 0}, 
 \    3   /    \    3   /    \    3   /    \    3   /            

  {t = Pi}, {t = 2 Pi}

 

restart
ode := 1/r*diff(r*u(r), r)=2*a:
dsolve(ode, u(r));
                                 2      
                              a r  + _C1
                       u(r) = ----------
                                  r     

You can also use expand(%) after the dsolve command if you prefer

2 real solutions and 4 complex ones :
 

restart:

e1 := x^2+2*m*x-6*x-m^3;
e2 := (x-r)*(x-r^2);
sol := solve([coeffs(expand(e1-e2), x)], [m, r])

-m^3+2*m*x+x^2-6*x

 

(x-r)*(-r^2+x)

 

[[m = 2, r = -2], [m = -3, r = 3], [m = -(1/2)*RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)^2-(1/2)*RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)+3, r = RootOf(_Z^4+4*_Z^3-5*_Z^2-24*_Z+36)]]

(1)

# check
map(s -> simplify(eval(e1-e2, s)), sol)

[0, 0, 0]

(2)

simplify([allvalues(sol[3])])

[[m = 1/4+((3/4)*I)*3^(1/2)+(1/4)*(25-(4*I)*3^(1/2))^(1/2)-((1/4)*I)*3^(1/2)*(25-(4*I)*3^(1/2))^(1/2), r = -1+((1/2)*I)*3^(1/2)+(1/2)*(25-(4*I)*3^(1/2))^(1/2)], [m = 1/4+((3/4)*I)*3^(1/2)-(1/4)*(25-(4*I)*3^(1/2))^(1/2)+((1/4)*I)*3^(1/2)*(25-(4*I)*3^(1/2))^(1/2), r = -1+((1/2)*I)*3^(1/2)-(1/2)*(25-(4*I)*3^(1/2))^(1/2)], [m = 1/4-((3/4)*I)*3^(1/2)-(1/4)*(25+(4*I)*3^(1/2))^(1/2)-((1/4)*I)*3^(1/2)*(25+(4*I)*3^(1/2))^(1/2), r = -1-((1/2)*I)*3^(1/2)-(1/2)*(25+(4*I)*3^(1/2))^(1/2)], [m = 1/4-((3/4)*I)*3^(1/2)+(1/4)*(25+(4*I)*3^(1/2))^(1/2)+((1/4)*I)*3^(1/2)*(25+(4*I)*3^(1/2))^(1/2), r = -1-((1/2)*I)*3^(1/2)+(1/2)*(25+(4*I)*3^(1/2))^(1/2)]]

(3)

# check
map(s -> simplify(eval(e1-e2, s)), %)

[0, 0, 0, 0]

(4)

 


 

Download AllSolutions.mw

restart
phi := (1+sqrt(5))/2;

_phi := `#mi("&phi;")`;


plots[pointplot]([seq([n, sin(n*phi)], n = 1000 .. 2000)], symbol = point, axes = boxed, labels = [n, sin(_phi*n)], labeldirections = [horizontal, vertical])


(obtained with Maple 2015.2)

Did you notice that sec(theta - (Pi/2) * (2/Pi *(theta + Pi/4))) = sqrt(2) ?

sec.mw


 

restart:
eq := x^6+x^4*(1.02-0.50e-2*y)+x^2*(-0.4950e-2*y+0.625e-5*y^2)+0.625e-5*y^2=0:

s := [solve(eq, x)]:
# check the number of solutions
numelems(s);
                               6

# define a function (sol) of y onto s
sol := unapply(s, y):

# Arrange the 6 plots of real and imaginay parts of each solution 
# into a 2x3 matrix 

M := Matrix(
      2, 3,
      (i, j) -> plot([(Re, Im)(sol(y)[3*(i-1)+j])], y=-10..10, color=[blue, red], legend=['Re', 'Im'], title=cat("Solution ", 3*(i-1)+j))
     ):

# Display M
DocumentTools:-Tabulate(M, width=60);

A_solution.mw

You will find in the first part (2D input mode) of the attached file something which runs without error.

But the results are very strange and I suspect something's wrong in your FDTD scheme.

So the last part of this same file (1D input mode) contains a much simpler way to obtain the 5 points discretization of the laplacian (matrix form) and a time marching based on a simple forward euler scheme (dt=0.1).
The Explore command (could have been animate instead) will help you to verify qualitatively that the diffusion process seemsto perform correctly.
The peak at x=y=0 is the result of he overlap of the BCs.at this point (given the additive way you constructed your BCs you have counted twice the boundary corner (0,0)).

Download 2d_heat_eq_RK4_mmcdara.mw

Adaptation to a RK4 time scheme is easy

If you just want the extrema ( @vv I didn't know about the term exprema) you can also use this
(this method is slower than the classical vv's method but, at least in this case, it works well and requires almost no mathematical background)

restart:
with(Student[Calculus1]):
A := x^2 + 400/x:
ExtremePoints( A );
                           [   (2/3)]
                           [2 5     ]

See also the attached file
AAA.mw

4 5 6 7 8 9 10 Last Page 6 of 30