tomleslie

It is impossible...

Answered: tomleslie 13876

February 20 2020

1 1

for Maple (or anyone else) to determine real and imaginary parts without knowing whther your "parameters", ie {C, Lm, Rpr, Rs, g, lpr, ls, w} are real or imaginary.

If I ask Maple to assume that all of these are real, then I get the result in the attached

restart; interface(showassumed = 0); assume(C, real, Lm, real, Rpr, real, Rs, real, g, real, lpr, real, ls, real, w, real)

Z := Rs+I*w*ls+I*w*Lm*(Rpr/g+I*w*lpr)/(I*w*Lm+Rpr/g+I*w*lpr)-I/(C*w)

Rs+I*w*ls+I*w*Lm*(Rpr/g+I*w*lpr)/(I*w*Lm+Rpr/g+I*w*lpr)-I/(C*w)

(1)

Rs+w^2*Lm^2*Rpr*g/(Lm^2*g^2*w^2+2*Lm*g^2*lpr*w^2+g^2*lpr^2*w^2+Rpr^2)

(2)

w*ls+w*Lm*(Lm*g^2*lpr*w^2+g^2*lpr^2*w^2+Rpr^2)/(Lm^2*g^2*w^2+2*Lm*g^2*lpr*w^2+g^2*lpr^2*w^2+Rpr^2)-1/(C*w)

(3)

(4)

Download ReIm.mw

Spaghetti code...

Answered: tomleslie 13876

February 19 2020

3 2

You do not seem to undertsand the difference between a function of three variables - eg F(i, j, x), and multiple indexed functions of two variables, eg F[i](j,x). In fact the worksheet you supplied switches from one interpretation to the other multiple times, making your intention almosy impossible to understand.

In the attached I have defined a function of three variables which appears to give the answers you expect in the comparison tests.

You should bear in mind that this invoolved a lot of guesswork on my part, so it comes with no guarantees, but at least it is logical and consistent. I suggest you test extensively

Trying to write a program for the representation of a function F with a three variables i, j, x, where i, j are non-negative integers and x is real.

I have tried with F(i,j,x), F[i](j,x), F[i,j](x), but did not get the desired result. My objective is to get F as a function of i, j, x.

Here is my code....

>

(1)

>

Astar := proc (s) options operator, arrow; lambda/(lambda+s) end proc; Astardas := proc (j, x) options operator, arrow; if j = 0 then Astar(x) else eval(diff(Astar(s), `$`(s, j)), s = x) end if end proc

proc (j, x) options operator, arrow; if j = 0 then Astar(x) else eval(diff(Astar(s), `$`(s, j)), s = x) end if end proc

(2)

>

  F:= proc( i::nonnegint, j::nonnegint, x::realcons, N::posint:=100)
            local lambda:= 1.6,
                  mu:= 2.0,
                  phi:= 0.7,
                  alpha:= 0.6,
                  Astar:= s -> lambda/(lambda + s),
                  Astardas:= (j, x)-> `if`( j = 0,
                                            Astar(x),
                                            eval
                                            ( diff
                                              ( Astar(s),
                                                s$j
                                              ),
                                              s = x
                                            )
                                          ),
                  psi:= i-> `if`( i=N,
                                  1,
                                  (1 - alpha)/alpha^(N - i)
                                );
                  if   i<N
                  then if   j=0
                       then return 0;
                       elif x=phi # shouldn't really test equality of floats!
                       then return -Astardas(j + 1, x)*psi(i - 1)/(j + 1);
                       else return piecewise( j = 0,
                                              (Astar(x)*psi(i - 1) - psi(i))/(phi - x),
                                              (Astardas(j, x)*psi(i - 1) + j*F(i, j - 1, x))/(phi - x)
                                            )
                       fi;
                  else if x=phi # shouldn't really test equality of floats!
                       then return -Astardas(j + 1, x)*(psi(N) + psi(N - 1))/(j + 1)
                       else return piecewise( j = 0,
                                              (Astar(x)*(psi(N) + psi(N - 1)) - psi(N))/(phi - x),
                                              (Astardas(j, x)*(psi(N) + psi(N - 1)) + j*F(N, j - 1, x))/(phi - x)
                                            )
                       fi;
                  fi;
    end proc:

>

#
# Some comparisons
#
  psi:= i-> `if`( i=N,
                  1,
                  (1 - alpha)/alpha^(N - i)
                ):
  for j to 10 do
      F(10, j, phi), -Astardas(j + 1, phi)*psi(10 - 1)/(j + 1);
  end do;
  for j to 10 do
      F(100, j, phi), -Astardas(j + 1, phi)*(psi(N)+ psi(N - 1))/(j + 1);
  end do

(3)

>

Download badFun.mw

Becuase the keyword table...

Answered: tomleslie 13876

February 19 2020

1 0

appears in your equation.

Eq1 := .9*(diff(f(eta), eta, eta, eta))+.1*We*(diff(f(eta), eta, eta))*(diff(f(eta), eta, eta, eta))-(diff(f(eta), eta))^2+f(eta)*(diff(f(eta), eta, eta))^2+0.6e-1-.5*(diff(f(eta), eta))+.1*(table([w = 1.1]))(eta)*(1+beta1*(table([w = 1.1]))(eta))+.1*phi(eta)*N*(1+.1[c]*phi(eta)) = 0;

I'm not sure what you intend this usage to do?????

What were you trying to achieve??

Parameters not assigned...

Answered: tomleslie 13876

February 19 2020

1 0

In the code you post, you have

E = -.1; N := 0.1e-1; n = .1;

The first and last of these are equations not assignments. SO values for 'E' and 'n' willnot be substitued in your ODE system. This means your system will contain two additional "unknowns" - for which you need a couple of additional boundary conditions. Hence the error message.

If I convert the above to assignments, then everything seems to work. See the attached.

>

  restart:
  blt:= 5:
  lambda:= .1: delta:= .1:
  gamma2:= .1:
  lambda2:= .1: lambda[T]:= .7: lambda[C] = .1:
  Pr:= 2:
  Nb:= 0.1e-1:
  Nt:= 0.1e-1:
  Sc:= 0.1e-1:
  M:= .5: E:= -.1:
  N:= 0.1e-1: n:= .1:
  sigma[1]:= .1:
  alpha[0]:= 0: omega:= 0:
  Eq1:= diff(f(eta), eta, eta, eta)+gamma2*((diff(f(eta), eta, eta))*(diff(f(eta), eta, eta))-f(eta)*(diff(f(eta), eta, eta, eta, eta)))-(1+lambda2)*((diff(f(eta), eta))*(diff(f(eta), eta))-f(eta)*(diff(f(eta), eta, eta)))-(1+lambda2)*M*sin(omega)^2*(diff(f(eta), eta))+(1+lambda2)(lambda*(1+lambda[T]*theta(eta))*theta(eta)+lambda*N*(1+lambda[C]*phi(eta))*phi(eta))*cos(alpha[0]) = 0;
  Eq2:= diff(theta(eta), eta, eta)+Pr*f(eta)*(diff(theta(eta), eta))+Pr*Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Pr*Nt*((diff(theta(eta), eta))*(diff(theta(eta), eta))) = 0;
  Eq3:= diff(phi(eta), eta, eta)+Sc*f(eta)*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta, eta))/Nb-Sc*sigma[1]*(1+delta*theta(eta))^n*exp(-E/(1+delta*theta(eta))) = 0;
  Eq3:= diff(phi(eta), eta, eta)+Sc*f(eta)*(diff(phi(eta), eta))+Nt*(diff(theta(eta), eta, eta))/Nb-Sc*sigma[1]*(1+delta*theta(eta))^n*exp(-E/(1+delta*theta(eta))) = 0;
bcs1:= f(0) = 0, (D(f))(0) = 1, (D(f))(blt) = 0, (D(D(f)))(blt) = 0, theta(0) = 1, theta(blt) = 0, Nb*(D(phi))(0)+Nt*(D(theta))(0) = 0, phi(blt) = 0;

diff(diff(diff(f(eta), eta), eta), eta)+.1*(diff(diff(f(eta), eta), eta))^2-.1*f(eta)*(diff(diff(diff(diff(f(eta), eta), eta), eta), eta))-1.1*(diff(f(eta), eta))^2+1.1*f(eta)*(diff(diff(f(eta), eta), eta))+1.1 = 0

diff(diff(theta(eta), eta), eta)+2*f(eta)*(diff(theta(eta), eta))+0.2e-1*(diff(theta(eta), eta))*(diff(phi(eta), eta))+0.2e-1*(diff(theta(eta), eta))^2 = 0

diff(diff(phi(eta), eta), eta)+0.1e-1*f(eta)*(diff(phi(eta), eta))+1.000000000*(diff(diff(theta(eta), eta), eta))-0.1e-2*(1+.1*theta(eta))^.1*exp(.1/(1+.1*theta(eta))) = 0

(1)

>

  L:= [.1, .5, 1.0, 5.0, 10.0]:
  for k to 5 do
      R:= dsolve( eval({Eq1, Eq2, Eq3, bcs1}, M = L[k]),
                  [f(eta), theta(eta), phi(eta)],
                  numeric,
                  method = bvp[midrich],
                  maxmesh = 4096,
                  abserr = 0.1e-1,
                  output = listprocedure
              );
      X1 || k:= rhs(R[3]);
      X2 || k:= rhs(R[4]);
      X3 || k:= rhs(R[5]);
      Y1 || k:= rhs(R[6]);
      Y2 || k:= -rhs(R[7]);
      Z1 || k:= rhs(R[8]);
      Z2 || k:= -rhs(R[9])
  end do:

>	print([(X1 \|\| (1 .. 5))(0)]); print([(X2 \|\| (1 .. 5))(0)]); print([(Y1 \|\| (1 .. 5))(0)]); print([(Y2 \|\| (1 .. 5))(0)]); print([(Z1 \|\| (1 .. 5))(0)]); print([(Z2 \|\| (1 .. 5))(0)]);

(2)

>	plot([X1 \|\| (1 .. 5)], 0 .. blt, labels = [eta, (D(f))(eta)], thickness = 1, color = black); plot([Y1 \|\| (1 .. 5)], 0 .. blt, labels = [eta, theta(eta)], thickness = 1, color = black); plot([Z1 \|\| (1 .. 5)], 0 .. blt, labels = [eta, phi(eta)], thickness = 1, color = black);

>

Download odeProb.mw

Or...

Answered: tomleslie 13876

February 17 2020

0 1

use

simplify(expr, power, symbolic)

use partition()...

Answered: tomleslie 13876

February 14 2020

1 1

as in the attached

Download cpart.mw

Just use diff()...

Answered: tomleslie 13876

February 14 2020

1 0

See the examples in the attached

>	restart; n:=2: M1:=Matrix( 2, 2, [ [ sin(t), t^2 ], [ ln(t), cos(t)] ] ); M2:=diff(M1, t$n);

$Matrix(2, 2, {(1, 1) = sin(t), (1, 2) = t^2, (2, 1) = ln(t), (2, 2) = cos(t)})$

Matrix(%id = 18446744074376952638)

(1)

>	restart; M1:=Matrix( 2, 2, [ [ sin(t), t^2 ], [ ln(t), cos(t)] ] ); M2:=diff(M1, t$n); eval(M2, n=2);

$Matrix(2, 2, {(1, 1) = sin(t), (1, 2) = t^2, (2, 1) = ln(t), (2, 2) = cos(t)})$

$diff(Matrix(2, 2, {(1, 1) = sin(t), (1, 2) = t^2, (2, 1) = ln(t), (2, 2) = cos(t)}), [`$`(t, n)])$

Matrix(%id = 18446744074376938910)

(2)

>

Download matDiff.mw

I'm not completely clear...

Answered: tomleslie 13876

February 13 2020

0 1

what the problem is here, but maybe one of the two options in the attached is what you are after?

Download getRoot.mw

Another alternative...

Answered: tomleslie 13876

February 12 2020

1 2

using Kitonum's (better) method for generating "grid" lines"

  display
  ( [ plot3d
      ( [ [x,t,0], [x,t,1], [0,x,t], [1,x,t], [x,0,t], [x,1,t] ],
          x=0..1,
          t=0..1,
          style=line,
          grid=[6,6],
          color=grey
      ),
      plot3d
      ( UU(x,t),
        x=0..xbas,
        t=0..tbas,
        color=cyan,
        transparency=0.75
      ),
      pointplot3d
      ( [ seq
          ( seq
            ( [ i, j, u_exact(i,j)],
              i=0..1.0,0.2
            ),
            j=0..1.,0.2
          )
        ],
        symbol=solidsphere,
        symbolsize=24,
        color=red
      )
    ]
  )

which produces

A lot depends...

Answered: tomleslie 13876

February 12 2020

1 2

on just how hard you are prepared to work to put all the pieces together in a single plot.

In addition to the solution given by Kitonum, see the final plot in the attached, which is pretty close to the second example in your original post

>

>	restart: with(orthopoly): with(LinearAlgebra): with(plots): interface(rtablesize=20):

>	k:=2: M:=2: teta:=1: beta:=1: K:=2^(k-1): N:=K*M:

>

lambda:=(m,teta,beta)->sqrt((2*m+teta+beta+1)*GAMMA(2*m+teta+beta+1)*m!/(2^(teta+beta+1)*GAMMA(m+teta+1)*GAMMA(m+beta+1)));
psii:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K, 2^(k/2) *lambda(m,teta,beta)*simplify(JacobiP(m,teta,beta,2^(k)*x-2*n+1)), 0);
local i,j: psi:=(x)->Array([seq(seq(psii(i,j,x),j=0..M-1),i=1..K)] ):
w:=(n,x)->(1+x)^(1/gama-1):

omega:=(n,x)->w(n,2^(k)*x-2*n+1):

lambda := proc (m, teta, beta) options operator, arrow; sqrt((2*m+teta+beta+1)*GAMMA(2*m+teta+beta+1)*factorial(m)/(2^(teta+beta+1)*GAMMA(m+teta+1)*GAMMA(m+beta+1))) end proc

proc (n, m, x) options operator, arrow; piecewise((n-1)/K <= x and x <= n/K, 2^((1/2)*k)*lambda(m, teta, beta)*simplify(JacobiP(m, teta, beta, 2^k*x-2*n+1)), 0) end proc

(1.1)

>	psi:=t->Matrix(N,1,[seq(seq(psii(i,j,t),j=0..M-1),i=1...K)] ):

>	for i from 1 to N do X[i]:=evalf((2i-1)/((2^k)M)): end do: for i from 1 to N do T[i]:=evalf((2i-1)/((2^k)M)): end do: for i from 1 to N do r[i]:=evalf(psi(T[i])): end do: Phi_mxm:=Matrix([seq(r[i],i=1...N)]);

Matrix(4, 4, {(1, 1) = 1.732050808, (1, 2) = 1.732050808, (1, 3) = 0., (1, 4) = 0., (2, 1) = -3.872983346, (2, 2) = 3.872983346, (2, 3) = 0., (2, 4) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.732050808, (3, 4) = 1.732050808, (4, 1) = 0., (4, 2) = 0., (4, 3) = -3.872983346, (4, 4) = 3.872983346})

(1.2)

>

P:=proc(k,M,nn) local PB,m,i,p,j,xi;
m:=M*(2^(k-1)):
xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
PB:=Matrix(m,m):
for i from 1 to m do
p:=0:
for j from 1 to m do
if i=j then PB[i,j]:=1;
fi:
if i<j then p:=p+1:
PB[i,j]:= xi(p,nn);
fi:
end do;
p:=1:
end do:
PB:=1/m^nn/GAMMA(nn+2)*PB;

return PB;
end proc:
PP:=alpha->Phi_mxm.P(k,M,alpha).MatrixInverse(Phi_mxm);

(1.3)

>	f1:=unapply(1/(1+exp(x))^2,x): g1:=unapply(1/(1+exp(-5t))^2,t); g2:=unapply(1/(1+exp(1-5t))^2,t); alpha:=1: #Analytical solution u_exact:=(x,t)-> 1/(1+exp(x-5*t))^2 :

>

U:=Matrix( N,N,symbol=u);
wwxx:= diff(f1(x),x,x)+(psi(x)^+.U.PP(alpha).psi(t))(1):
wwt:= diff(g1(t),t)+x*(diff(g2(t),t)-diff(g1(t),t)-(psi(1)^+.PP(2)^+.U.psi(t))+psi(x)^+.PP(2)^+.U.psi(t))(1) :
ww:= f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1) :
MC:=unapply(wwt -wwxx -6*ww .(1-ww ),x,t):

g1 := proc (t) options operator, arrow; 1/(1+exp(-5*t))^2 end proc

g2 := proc (t) options operator, arrow; 1/(1+exp(1-5*t))^2 end proc

Matrix(%id = 18446744074365962950)

(1)

>	PDEson:=Matrix(N,N): for i from 1 to N do for j from 1 to N do PDEson(i,j):=simplify( evalf( MC(X[i],T[j]) ) )=0: end do: end do:

>	coz:=fsolve({seq(seq(PDEson(i,j),i=1..N ),j=1..N )}): U:=subs(op(coz),U):

>	#Numerical Solution UU:=(x,t)->evalf(f1(x)+g1(t)-g1(0)+x*(g2(t)-g2(0)-g1(t)+g1(0)-(psi(1)^+.PP(2)^+.U.PP(alpha).psi(t)))(1)+(psi(x)^+.PP(2)^+.U.PP(alpha).psi(t))(1)):

>

U_ERR:=(x,t)->abs(UU(x,t)-u_exact(x,t)):
xbas:=1:
tbas:=1:
with(plots):
#printf("Numerical Solution");
g11:=plot3d(UU(x,t),x=0..xbas,t=0..tbas,color=gray,title='NumericalSolution');
#printf("Exact Solution");
g22:=plot3d(u_exact(x,t),x=0..xbas,t=0..tbas,title='ExactSolution');
#printf("Both exact and numerical solution");
display(g11,g22, title='Both_exact_and_numerical_solution');
#printf("Error");
plot3d(U_ERR(x,t),x=0..xbas,t=0..tbas,title='Error');

>

  with(plottools):
  display
  ( [ seq
      ( seq
        ( line
          ( [j,i,0],
            [j,i,1],
            color=grey
          ),
          i=0..1,0.2
        ),
        j=0..1
      ),
      seq
      ( seq
        ( line
          ( [i,j,0],
            [i,j,1],
            color=grey
          ),
          i=0..1,0.2
        ),
        j=0..1
      ),
      seq
      ( seq
        ( line
          ( [j,0,i],
            [j,1,i],
            color=grey
          ),
          i=0..1,0.2
        ),
        j=0..1
      ),
      seq
      ( seq
        ( line
          ( [0,j,i],
            [1,j,i],
            color=grey
          ),
          i=0..1,0.2
        ),
        j=0..1
      ),
      plot3d
      ( UU(x,t),
        x=0..xbas,
        t=0..tbas,
        color=gray
      ),
      pointplot3d
      ( [ seq
          ( seq
            ( [ i, j, u_exact(i,j)],
              i=0..1.0,0.2
            ),
            j=0..1.,0.2
          )
        ],
        symbol=solidsphere,
        symbolsize=24,
        color=red
      )
    ]
  );

>

Download 3DGR.mw

or cat()...

Answered: tomleslie 13876

February 11 2020

0 0

Your fundamental problem is that in your original code the character 'i' appears within a string. As a result it will always be a simple string character 'i' with no reference to the loop index. After all, would you expect the character 'i' in the string "Zubair" to be replaced by the integer value of the loop variable??

The conventional way to handle this is to use the cat() command (I recommend reading the help page) as in the attached

Download catuse.mw

Identifies the variable...

Answered: tomleslie 13876

February 11 2020

1 1

The subscript on the 'D' simply identifies (by position) the variable with respect to which differentiation is being performed, so

D₁(f)(x,y) is diff( f(x,y), x), and

D₂(f)(x,y) is diff( f(x,y), y)

see the attached for details

>

#
# Start with a "trivial" example to show
# that D-subscript-1 just means differentiation
# with respect to the *first* argument of the
# function and D-subscript-2 means
# differentiation with respect to the second
# argument of the function
#
  diff(F(x,y), x);
  convert(%, D);
  diff(F(x,y), y);
  convert(%, D);

(1)

>	# # OP's example combines the above notation with a # "chain-rule" example, because # # diff( F(x, g(y)), y) is, (by the chain rule) # # diff( F(), g)* diff(g(), y) # # which can be seen from # diff(F(x, g(y)), y); convert(%, diff);

(eval(diff(F(x, t1), t1), {t1 = g(y)}))*(diff(g(y), y))

(2)

>

Download diffDef.mw

An alternative...

Answered: tomleslie 13876

February 10 2020

0 0

Don't "unlist" your list (ie just use 'd' rather than 'd[]') before printing with %a, as in the attached

Download pList.mw

Not completely clear...

Answered: tomleslie 13876

February 10 2020

0 0

(to me at least) which curves you actually want (or the form in which you want them).

Maybe one of the two options in the attached?

>

  restart:
  with(plots):
  p := 1: p1 := -2:
  para := -2*p*x+y^2 = 0:
  para1 := -2*p1*x+y^2 = 0:
  t := y-m*x-p/(2*m) = 0:
  t1 := y+x/m+(1/2)*p1*m = 0:
  sol := op(solve([t, t1], [x, y]));
  plot( [ [rhs~(sol)[], m=0..10],
          [rhs~(sol)[], m=-10..0]
        ],
        labels=["xval", "yval"]
      );
  spacecurve( [m, rhs(sol[1]), rhs(sol[2])],
              m=-10..10,
              thickness=4,
              color=red,
              labels=["m", "xval", "yval"]
            );

[x = (1/2)*(2*m^2-1)/(m^2+1), y = (1/2)*(2*m^4+1)/(m*(m^2+1))]

>

Download pltver.mw

Hmmmm...

Answered: tomleslie 13876

February 10 2020

1 0

Some observations (and a pretty crude worksheet), because I'm pressed for time at the moment

OP generates an equation f( w, y(w), m, _C1 )=0 which arises as an implicit solution of an ODE. 'm' is apparently a 'parameter', and _C1 is an integration constant
OP produces a new equation by setting [m=1, _C1=0]. This equation is of the form g( w, y(w) )=0, but it is not possible to isolate the dependent function 'y(w)'.
Under certain circumstances it would be possible to plot the solution of the equation obtained in (2) above by using the implicitplot() command. However a little experimentation shows that for 'w' real, 'y(w)' is always(?) complex - and implicitplot() does't like comple solutions!
An alternative is to use fsolve(..., complex) to solve the equation generated in (2) above, producing complex-valued 'y(w)' for any given value of 'w'.
One can then use complexplot() to plot these solutions - although note that
1. This plots Im( y(w) ) against Re( y(w) ) with 'w' being a parameter for the plot
2. I could only make this work for negative values of the variable 'w'
It is also possible to use spacecurve() to plot the curve with the components [w, Re(y(w)), Im(y(w))] - this takes a while

See the attached

>

diff(y(w), w)+(2*3^(1/2)*Pi(y(w)^2+m^2*w^2)^(1/2)*y(w)+m^2*w)/y(w) = 0

(1)

>

$[{ln(w)+(1/4)*ln(-m^4-2*m^2*y(w)^2/w^2-y(w)^4/w^4+12*y(w)^2*Pi/w^2)-(3/2)*arctan((1/4)*(-2*m^2-2*y(w)^2/w^2+12*Pi)/(3*Pi*m^2-9*Pi^2)^(1/2))*Pi/(3*Pi*m^2-9*Pi^2)^(1/2)+(1/4)*ln(2*Pi^(1/2)*3^(1/2)*y(w)/w-y(w)^2/w^2-m^2)-(3/2)*arctan((1/2)*(2*Pi^(1/2)*3^(1/2)-2*y(w)/w)/(m^2-3*Pi)^(1/2))/((3*m^2-9*Pi)/Pi)^(1/2)-(1/4)*ln(2*Pi^(1/2)*3^(1/2)*y(w)/w+m^2+y(w)^2/w^2)+(3/2)*arctan((1/2)*(2*y(w)/w+2*Pi^(1/2)*3^(1/2))/(m^2-3*Pi)^(1/2))/((3*m^2-9*Pi)/Pi)^(1/2)-_C1 = 0}]$

(2)

>

#
# Solve 'Ans' for y(w) in terms of w. The solution is
# always(?) complex. Furthermore 'fsolve' only returns
# a solution for w<0. A few examples below.
#
  eq2:=eval( Ans[1,1],
             [ _C1=0, m=1, y(w)=Y ]
           ):
  fsolve( eval(eq2, w=-2 ), Y, complex);
  fsolve( eval(eq2, w=-1 ), Y, complex);
  fsolve( eval(eq2, w=0 ), Y, complex);
  fsolve( eval(eq2, w=1 ), Y, complex);
  fsolve( eval(eq2, w=2 ), Y, complex);
#
# (Complex)plot y(w) as a function of w over the
# range w=-2.0..-0.01
#
# NB this takes a couple of minutes
#
  plots:-complexplot
         ( 'fsolve'
           ( eval(eq2, w=zz),
             Y,
             complex
           ),
           zz=-2..-0.01
         );

Error, (in ln) numeric exception: division by zero

fsolve((1/4)*ln(-Y^4+12*Pi*Y^2-2*Y^2-1)-(3/2)*arctan((1/4)*(-2*Y^2+12*Pi-2)/(-9*Pi^2+3*Pi)^(1/2))*Pi/(-9*Pi^2+3*Pi)^(1/2)+(1/4)*ln(2*Pi^(1/2)*3^(1/2)*Y-Y^2-1)-(3/2)*arctan((1/2)*(2*Pi^(1/2)*3^(1/2)-2*Y)/(1-3*Pi)^(1/2))/((-9*Pi+3)/Pi)^(1/2)-(1/4)*ln(2*Pi^(1/2)*3^(1/2)*Y+1+Y^2)+(3/2)*arctan((1/2)*(2*Y+2*Pi^(1/2)*3^(1/2))/(1-3*Pi)^(1/2))/((-9*Pi+3)/Pi)^(1/2) = 0, Y, complex)

fsolve(ln(2)+(1/4)*ln(-1-(1/2)*Y^2-(1/16)*Y^4+3*Y^2*Pi)-(3/2)*arctan((1/4)*(-2-(1/2)*Y^2+12*Pi)/(-9*Pi^2+3*Pi)^(1/2))*Pi/(-9*Pi^2+3*Pi)^(1/2)+(1/4)*ln(Pi^(1/2)*3^(1/2)*Y-(1/4)*Y^2-1)-(3/2)*arctan((1/2)*(-Y+2*Pi^(1/2)*3^(1/2))/(1-3*Pi)^(1/2))/((-9*Pi+3)/Pi)^(1/2)-(1/4)*ln(Pi^(1/2)*3^(1/2)*Y+1+(1/4)*Y^2)+(3/2)*arctan((1/2)*(2*Pi^(1/2)*3^(1/2)+Y)/(1-3*Pi)^(1/2))/((-9*Pi+3)/Pi)^(1/2) = 0, Y, complex)

>

#
# Using a similar approach, use spacecurve to produce
# a plot of the curve [w, Re(y(w)), Im(y(w))
#
# NB this takes a couple of minutes
#
  plots:-spacecurve
         ( [ zz,
             Re( 'fsolve'
                 ( eval(eq2, w=zz),
                   Y,
                   complex
                 )
               ),
             Im( 'fsolve'
                 ( eval(eq2, w=zz),
                   Y,
                   complex
                 )
                )
           ],
           zz=-2..-0.01,
           labels=[ "w", "Re(y(w))", "Im(y(w))"],
           thickness=4,
           color=red
         );

>

Download compPlot.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

13876 Reputation

20 Badges

MaplePrimes Activity

These are answers submitted by tomleslie

It is impossible...

Spaghetti code...

Becuase the keyword table...

Parameters not assigned...

Or...

use partition()...

Just use diff()...

I'm not completely clear...

Another alternative...

A lot depends...

or cat()...

Identifies the variable...

An alternative...

Not completely clear...

Hmmmm...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again