acer

32303 Reputation

29 Badges

19 years, 310 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You can wrap the result from seq in square brackets, in order to make a list from the sequence of numbers.

  L := [seq(a(n), n = 1 .. 100)];
  plots:-listplot(L);

 

Do you mean like this?

seq([seq([a[i,j],b[(j-1)*3+i]],i=1..3)], j=1..3);

   [[a[1, 1], b[1]], [a[2, 1], b[2]], [a[3, 1], b[3]]],

   [[a[1, 2], b[4]], [a[2, 2], b[5]], [a[3, 2], b[6]]],

   [[a[1, 3], b[7]], [a[2, 3], b[8]], [a[3, 3], b[9]]]

Or, if with slightly different bracketing,

seq(seq([a[i,j],b[(j-1)*3+i]],i=1..3), j=1..3);

   [a[1, 1], b[1]], [a[2, 1], b[2]], [a[3, 1], b[3]],
   [a[1, 2], b[4]], [a[2, 2], b[5]], [a[3, 2], b[6]],
   [a[1, 3], b[7]], [a[2, 3], b[8]], [a[3, 3], b[9]]

This runs in Maple 18.02 (since your Questions are usually restricted to that old version).

restart;

kernelopts(version);

`Maple 18.02, X86 64 LINUX, Oct 20 2014, Build ID 991181`

H := a__1*exp(I*(-Omega*t+k*x+z*k[1]))
     +a__2*exp(-I*(-Omega*t+k*x+z*k[1]));

a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1]))

expr := I*(diff(H, z))+diff(H, x, x)+diff(H, t, t)
        +R*(H+conjugate(H))+R^2*(H+conjugate(H))*H;

I*(I*a__1*k[1]*exp(I*(-Omega*t+k*x+z*k[1]))-I*a__2*k[1]*exp(-I*(-Omega*t+k*x+z*k[1])))-a__1*k^2*exp(I*(-Omega*t+k*x+z*k[1]))-a__2*k^2*exp(-I*(-Omega*t+k*x+z*k[1]))-a__1*Omega^2*exp(I*(-Omega*t+k*x+z*k[1]))-a__2*Omega^2*exp(-I*(-Omega*t+k*x+z*k[1]))+R*(a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1]))+conjugate(a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1]))))+R^2*(a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1]))+conjugate(a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1]))))*(a__1*exp(I*(-Omega*t+k*x+z*k[1]))+a__2*exp(-I*(-Omega*t+k*x+z*k[1])))

expr2 := subs(-Omega*t+k*x+z*k[1]=W(x),expr):

collect(simplify(convert(evalc(expr2),exp)),exp,u->simplify(u,size)):

subs(W(x)=-Omega*t+k*x+z*k[1],subsindets(%,specfunc(exp),expand)):
ans:=subsindets(%,And(specfunc(exp)^(integer),
                       satisfies(u->op(2,u)<0)),
                u->exp(-op([1,..],u))^abs(op(2,u)));

R^2*a__2*(a__1+a__2)*(exp(-I*(-Omega*t+k*x+z*k[1])))^2+((-Omega^2-k^2+R+k[1])*a__2+R*a__1)*exp(-I*(-Omega*t+k*x+z*k[1]))+((-Omega^2-k^2+R-k[1])*a__1+R*a__2)*exp(I*(-Omega*t+k*x+z*k[1]))+R^2*a__1*(a__1+a__2)*(exp(I*(-Omega*t+k*x+z*k[1])))^2+R^2*(a__1+a__2)^2

simplify(evalc(ans-expr));

0

Download collect_exp2_18.mw

In the above code the term
   exp(-I*(-Omega*t+k*x+z*k[1]))
is deliberately kept as you gave it, ie. not turned into
   exp(I*(Omega*t-k*x-z*k[1]))
or
   1/exp(I*(-Omega*t+k*x+z*k[1]))

And similarly for the other exp call.

Also, a term like
   exp(-I*(-Omega*t+k*x+z*k[1]))^2
is deliberately kept as is, and not turned into
   exp(-2*I*(-Omega*t+k*x+z*k[1]))

You have not informed us as to which variables are to be considered as purely real. So I am guessing that all the indeterminates can be taken as purely real. You originally omitted similar information in another of your recent Questions. This is not helpful.

You asked why your loop approach cannot be replaced with seq.

So I will show you how your seq attempt can be easily adjusted to work, using the prefix form of `+=`.

Using your N, V, and sumidx,

roll:=rand(1..9):
N := sort([seq(seq(10*(2+jdx)+roll(),jdx=1..5),idx=1..10)]):
V := 10*seq(roll()+10*trunc(N[idx]/10),idx=1..nops(N)):
sumidx := [seq(floor(N[idx]/10)-2,idx=1..nops(N))];

You can now proceed with,

S := 0 *~convert(convert(trunc~(N/10),set),list):
seq(`+=`(S[ sumidx[idx] ],V[idx]), idx=1..nops(N)):
S;

     [3420, 4570, 5480, 6620, 7320]

which agrees with your loop approach,

S := 0 *~convert(convert(trunc~(N/10),set),list):
for idx to nops(sumidx) do
    S[ sumidx[idx] ] += V[idx]:
end do:
S;

     [3420, 4570, 5480, 6620, 7320]

A slightly wider pair of ranges for the parameters (along with a restricted vertical view) can help avoid the impression that the two surfaces are equal along a certain direction.

You can also play around with the stylistic effects.

plot3d([exp(t^2+10*u-1), -121+22*t+10*u, [[11, -12, 1]]],
       t=8..12, u=-15..-6, view=-50..50,
       color=[grey,green,red], style=[patch,surface,point],
       symbol=solidcircle, symbolsize=10,
       lightmodel=Light2, transparency=[0.4,0.0,0.0]);

It can get tricky finding a decent view, since a wider range (ie. larger t) alongside the restricted view can cause the GUI's plot render to fail to show one of the surfaces.

Here is one way to avoid calls to TF from happening (prematurely), ie. unless its arguments are not all of type realcons.

restart;

TF := proc(x, x0, L, k, alpha_0, alpha_L)
  if not [args]::list(realcons) then
    return 'procname'(args);
  end if;
  if x0 <= x then
    (abs(cos(k*(x - L) - alpha_L)
    *cos(k*x0 + alpha_0)/(cos(alpha_L)*cos(k*x0 + alpha_0))));
  else
    (abs(cos(k*(x0 - L) - alpha_L)
    *cos(k*x + alpha_0)/(cos(alpha_L)*cos(k*x0 + alpha_0))));
  end if;
end proc:

Optimization:-Minimize(TF(x, 0, 0.03, 55.11566060, Pi/2, Pi/4), x = 0 .. 0.03);

[HFloat(4.281548719743726e-7), [x = HFloat(0.015749994509751696)]]

plot(TF(x, 0, 0.03, 55.11566060, Pi/2, Pi/4), x = 0 .. 0.03, size=[500,300]);

Download Min_ex.mw

Here are some ideas, which you should be able to adjust.

I find that for animations of compound plots it's useful to create a procedure (below, F) which can be used to quickly generate and test each frame -- calling it with only a value for the animating parameter. Eg,  F(0.5) etc.

Then you can more quickly adjust that procedure, as needed.

And then the call to plots:-animate is very simple as a last step.

restart;

f1 := x->a*x^2:
f2 := x->a*x^2 + 0.3:

expr1 := f1(x);
expr2 := f2(x);

a*x^2

a*x^2+.3

F := proc(A) local x1,x2,y1,y2;
  uses plots:
  x1,x2 := 1.0, A/2;
  y1 := eval(expr1,[a=A, x=x1]);
  y2 := eval(expr2,[a=A, x=x2]);
  display(
    plot(eval([expr1, expr2],a=A), x=-1 .. 1,
         'color'=["Green","Red"], 'thickness'=3),
    pointplot([[x1,y1], [x2,y2]],
              'symbolsize'=15, 'color'=["Green","Red"]),
    textplot([[x1, y1, sprintf("f1(1.0)=%.4f", y1)],
              [x2, y2, sprintf("f2(a/2)=%.4f", y2)]],
             'align'=[':-below',':-right']));
end proc:

 

plots:-animate(F, [a], a = 0.2 .. 1.2);

 

Download plotting_anim_textplot.mw

You are accidentally passing a sequence of scalar expressions to the simplify command, since that is what dsolve is returning.

You could correct that by wrapping those results from dsolve in a list, eg.

   simplify([dsolve(DE, explicit)]);

When you pass such a sequence of arguments to simplify then some of the arguments are taken as being options to the simplify command. But your goal appears to be simplification of all that dsolve returns (which could be a sequence of multiple solutions). So you could pass them all within a list or set.

Here's one way, using your operators assigned to f1 and f2.

f1:=x->a*x^2:
f2:=x->a*x^2+1:

Explore(plot(eval([f1(x),f2(x)],a=A),
             x=-1 .. 1, y = -3 .. 3),
        parameters = [[A = -1.0 .. 1.0, label=a]]);

It'd be easier (as Kitonum has already shown) if you simply used expressions for f1 and f2, instead of operators. That also works with assigned names f1 and f2.

f1:=a*x^2:
f2:=a*x^2+1:

Explore(plot([f1,f2],
             x=-1 .. 1, y = -3 .. 3),
        a = -1.0 .. 1.0);

Are x, t, and lambda known to be purely real?

It might be reasonably taken that the other names lambda1, epsilon1, and epsilon2 are not purely real, because they were entered stand-alone in calls to conjugate. It seems weird that you'd enter those if you didn't intend it so. If this is true then using only evalc (without further assumptions on these other names) is not appropriate.

The following also works in Maple 18.02.

restart;
alias(phi = phi(x, t), chi = chi(x, t), psi = psi(x, t), rho = rho(x, t)):

A := -phi*(2*epsilon1*conjugate(epsilon1)*psi*conjugate(epsilon2)*conjugate(phi)*epsilon2*conjugate(psi)+epsilon1^2*conjugate(epsilon1)^2*psi*conjugate(phi)*conjugate(psi)+chi*conjugate(epsilon2)*epsilon2*conjugate(psi)*conjugate(rho)+epsilon1*rho*conjugate(epsilon1)*conjugate(phi)*conjugate(rho)+epsilon2^2*psi*conjugate(epsilon2)^2*conjugate(phi)*conjugate(psi))*lambda1-conjugate(lambda1)*conjugate(chi)*(epsilon2*psi*rho*conjugate(epsilon2)*conjugate(phi)+chi*conjugate(epsilon1)*psi*epsilon1*conjugate(psi)+chi*rho*conjugate(rho)):

expr := subs({chi = exp((I*lambda*(1/2))*x-I*t/lambda),
              phi = exp(-(I*lambda*(1/2))*x+I*t/lambda),
              psi = exp(-(I*lambda*(1/2))*x-I*t/lambda),
              rho = exp((I*lambda*(1/2))*x+I*t/lambda)}, A):

simplify(expand(expr), size) assuming lambda::real, t::real, x::real;

    -(epsilon1*conjugate(epsilon1) + epsilon2*conjugate(epsilon2) + 1)
    *(lambda1*epsilon1*conjugate(epsilon1)
      + lambda1*conjugate(epsilon2)*epsilon2
      + conjugate(lambda1));

Slightly shorter and lighter is,

   factor(expand(expr)) assuming lambda::real, t::real, x::real;

You could also follow that up with,

   map(simplify,%);

to get the products of the form name*conjugate(name) turned into abs(name)^2. Or you could issue convert(%,abs) as a last step using the code above.

restart;

A := Vector(3,symbol=v);

Vector(3, {(1) = v[1], (2) = v[2], (3) = v[3]})

v[1] := 5;

5

Access to A[1] gets the evaluation, and returns the
running value of v[1].

A[1];

5

eval~(A);

Vector(3, {(1) = 5, (2) = v[2], (3) = v[3]})

rtable_eval(A);

Vector(3, {(1) = 5, (2) = v[2], (3) = v[3]})

rtable_eval(A, inplace);

Vector(3, {(1) = 5, (2) = v[2], (3) = v[3]})

A;

Vector(3, {(1) = 5, (2) = v[2], (3) = v[3]})

eqns := { 2 = c[2] + a[1],
          2 = 2*c[3] + a[2],
         -2 = -3*c[1] + a[3],
         -1 = -2*c[2] + a[1],
         -1 = -c[3] + a[2],
         1 = a[3]}:

res := solve(eqns);

{a[1] = 1, a[2] = 0, a[3] = 1, c[1] = 1, c[2] = 1, c[3] = 1}

V1 := Vector(3,symbol=a):
V2 := Vector(3,symbol=c):

eval(V1, res), eval(V2, res);

Vector(3, {(1) = 1, (2) = 0, (3) = 1}), Vector(3, {(1) = 1, (2) = 1, (3) = 1})

assign(res);

V[1];

V[1]

eval~(V1), eval~(V2);

Vector(3, {(1) = 1, (2) = 0, (3) = 1}), Vector(3, {(1) = 1, (2) = 1, (3) = 1})

rtable_eval(V1), rtable_eval(V2);

Vector(3, {(1) = 1, (2) = 0, (3) = 1}), Vector(3, {(1) = 1, (2) = 1, (3) = 1})

 

Download vector_eval.mw

One alternative is to try and compute it numerically.

I get PathInt to return the inert integral directly, then try to compute that numerically. I do that in two ways. The first way just splits into real and imaginary components, and loosens the target accuracy.

restart;

with(VectorCalculus,PathInt):

ig := PathInt(4*x*y*(y^4+2*x*y-2)/sqrt((1-(2*x*y)^2+(-3*y^2+1)^2)*(1+(y^2-1)^2)),
              [x, y] = Path(`<,>`(-(1+sqrt(2))*cos(t)-(sqrt(2)-1)*sin(t),
                                  cos(t)-sin(t)), t = 0 .. 2*Pi), inert);

Int(4*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))*(cos(t)-sin(t))*((cos(t)-sin(t))^4+2*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))*(cos(t)-sin(t))-2)*((sin(t)*2^(1/2)-cos(t)*2^(1/2)+sin(t)+cos(t))^2+(-sin(t)-cos(t))^2)^(1/2)/((1-4*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))^2*(cos(t)-sin(t))^2+(-3*(cos(t)-sin(t))^2+1)^2)*(1+((cos(t)-sin(t))^2-1)^2))^(1/2), t = 0 .. 2*Pi)

normal( eval(op(1,ig),t=t+Pi) - op(1,ig) );

0

2 * evalf[15](Int(Re(op(1,ig)),t=0 .. Pi, epsilon=1e-6)
              + I*Int(Im(op(1,ig)),t=0 .. Pi, epsilon=1e-6)):
evalf[6](%);

10.2837-35.6961*I

Download some_path_int0.mw

The second way does it in steps, manually splitting the domain. I did it partly because I got once got some conflicting results when trying more straightforwardly, and I wanted to double-check.

I split the domain into piecewise continuous parts, then integrate the real and complex parts of the integrand separately (with a forced choice of numeric method) over each of those parts. I utilize the periodicity to keep t small (in the first Pi period) which seems to help numerically. Once the domain is split, this approach seems to allow quick and successful computation with a finer accuracy tolerance than I was able to take in my first approach.

restart;

ig := VectorCalculus:-PathInt(4*x*y*(y^4+2*x*y-2)/sqrt((1-(2*x*y)^2+(-3*y^2+1)^2)*(1+(y^2-1)^2)),
                              [x, y] = Path(`<,>`(-(1+sqrt(2))*cos(t)-(sqrt(2)-1)*sin(t),
                                                  cos(t)-sin(t)), t = 0 .. 2*Pi), inert);

Int(4*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))*(cos(t)-sin(t))*((cos(t)-sin(t))^4+2*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))*(cos(t)-sin(t))-2)*((sin(t)*2^(1/2)-cos(t)*2^(1/2)+sin(t)+cos(t))^2+(-sin(t)-cos(t))^2)^(1/2)/((1-4*(-(1+2^(1/2))*cos(t)-(2^(1/2)-1)*sin(t))^2*(cos(t)-sin(t))^2+(-3*(cos(t)-sin(t))^2+1)^2)*(1+((cos(t)-sin(t))^2-1)^2))^(1/2), t = 0 .. 2*Pi)

normal( eval(op(1,ig),t=t+Pi) - op(1,ig) );

0

dd := discont(op(1,ig),t):

dptsPi := evalf[25]([0, select(u->u::realcons and is(u>=0 and u<=Pi),
                   evalf[25](map[2](op@eval,dd,_Z1=~[0,1])))[],
                   Pi]);

[0., .5855649839981530717439709, 1.051491027553841209737523, 1.386655922771385485327700, 2.543319910635242826584368, 3.141592653589793238462643]

plots:-display(
  plot([Re,Im](op(1,ig)), t=0..2*Pi,
               color=[red,green], thickness=3),
  plots:-pointplot(`[]`~(dptsPi,0),color=blue,
                   symbol=solidcircle,symbolsize=15),
  size=[500,200],
  view=-30..30, xtickmarks=decimalticks);

fre := unapply(Re(op(1,ig)),t):
fim := unapply(Im(op(1,ig)),t):

opts := epsilon=1e-10, method=_d01ajc, maxintervals=10000:
2 * add( evalf[15](Int(fre, dptsPi[i]..dptsPi[i+1], opts)
                   + I*Int(fim, dptsPi[i]..dptsPi[i+1], opts)),
         i=1..nops(dptsPi)-1 ):
evalf[10](%);

10.28375512-35.69630236*I

 

Download some_path_int.mw

With luck I didn't mess it all up.

The only slow bit was the call to discont. (I didn't try fiddling with fdiscont.) But that only has to be determined once, after which that domain splitting and the periodicity allow arbitrary integration ranges to be handled quickly.

Here are two ways.

The goal itself seems dubious. If you want the body of procedure b to not depend on the changing nature of procedure a then why construct b in that that way? (Also, why do you have to use procedures instead of expressions, at all?)

restart;

a := x -> x + 1;

proc (x) options operator, arrow; x+1 end proc

a := subsop(3=[{op(3,eval(a)),':-inline'}[]][],
            eval(a));

proc (x) options inline, operator, arrow; x+1 end proc

b := x -> a(x) + x^2;

proc (x) options operator, arrow; x+1+x^2 end proc

a := x -> x + 5;

proc (x) options operator, arrow; x+5 end proc

b(x);

x^2+x+1

restart;

a := x -> x + 1;

proc (x) options operator, arrow; x+1 end proc

b := unapply(a(x) + x^2, x);

proc (x) options operator, arrow; x^2+x+1 end proc

a := x -> x + 5;

proc (x) options operator, arrow; x+5 end proc

b(x);

x^2+x+1

Download proc_s.mw

If you have total control over the initial construction of procedure a then you could also do the first variant as follows:

restart;

a := proc(x) option operator,arrow,inline;
       x + 1;
     end proc;

proc (x) options inline, operator, arrow; x+1 end proc

b := x -> a(x) + x^2;

proc (x) options operator, arrow; x+1+x^2 end proc

a := x -> x + 5;

proc (x) options operator, arrow; x+5 end proc

b(x);

x^2+x+1

Download proc_inline.mw

It's possible that this is only a toy example, and that you may have some more involved example for which you cannot sucessfully utilize the unapply approach because the nature of procedure b doesn't allow you to properly call it with an unassigned name like x as its argument.

If you have more restrictions, etc, then it would be better to show us.

If the procedures a and b are completely constructed by you (and not the return values of some computations) then using expressions instead of procedures seems so much more straightforward and usual.

Yes, that's one of the main points about "loading" a package using the with command. It rebinds the names of its exports.

You can still access the original "global" name by using the colon-dash syntax in your own entered/parsed code. For example,

restart;

f := x^11 + 2*x^9 + 2*x^8 + x^6 + x^5 + 2*x^3 + 2*x^2 + 1;

x^11+2*x^9+2*x^8+x^6+x^5+2*x^3+2*x^2+1

g := 2*x^10 + x^7 + 2*x^4 + x;

2*x^10+x^7+2*x^4+x

Gcd(f, g) mod 3;

x^9+2*x^6+x^3+2

with(Algebraic);

[Content, ConvertRootOf, Degree, Divide, Expand, ExtendedEuclideanAlgorithm, Gcd, Gcdex, GetAlgebraics, GreatestCommonDivisor, MakeMonic, Normal, PrimitivePart, PseudoDivision, Quotient, Reduce, Remainder, Resultant, Squarefree]

Gcd(f, g) mod 3;

x^6+1

:-Gcd(f, g) mod 3;

x^9+2*x^6+x^3+2

Download pkg_rebinding.mw

Note that instances of the name Gcd in procedures already (before the rebinding) saved in a .mla archives are not affected -- they were saved with the original name, whose value you have not reassigned. So merely loading the package doesn't change the value of the instances of the these global names inside the procedures of, say, stock Maple commands.

ps. I changed your query from a Post into a Question.

Is this the kind of this you're after?

I ran this in Maple 16.02, but perhaps you are running Maple 12?

The 3D plot seems to show that (for any value of Zstub in its given range) the maximal value of Gamma is attained at one of the end-points: either f=f1 or f=f2. And you seem to want to find the Zproc that minimizes that maximum. That seems to occur at Zstub=Z[N]. But, even so, you may still wish to compute that via the nested optimization, as opposed to merely observing where it occurs on the 3D plot.

restart:
with(plots):
Digits := 15:
with(linalg):
j := sqrt(-1):

RL := 26; #dB
Gamma_m := evalf(10^(-RL/20)); #unitless, min req'd
f1 := 1.9;
f2 := 6.3;
f0 := evalf((f1+f2)/2);
BW := evalf((f2-f1)/f0);
ZS := 50/16;
ZL := 50;

26

0.501187233627272e-1

1.9

6.3

4.10000000000000

1.07317073170732

25/8

50

Chebyshev

 

# A NEW WAY OF FINDING 'N'
theta_m := 'theta_m':
N := 'N':
theta_m := solve(BW=2-4*theta_m/Pi,theta_m);
BW := evalf(2-4*theta_m/Pi);  # used at the end to find shoulders (could maybe use 0..f0 and f0..f2+f0/2 instead
N := solve(sec(theta_m)=cosh((1/N)*arccosh(abs(ln(ZL/ZS)/(2*Gamma_m)))),N);
printf("N = %g yields RL = %g\n", floor(N),
  20*log10(abs(fsolve(sec(theta_m)=cosh((1/floor(N))*arccosh(abs(ln(ZL/ZS)/(2*Gamma)))),Gamma=Gamma_m))));
printf("N = %g yields RL = %g\n", ceil(N),
  20*log10(abs(fsolve(sec(theta_m)=cosh((1/ceil(N))*arccosh(abs(ln(ZL/ZS)/(2*Gamma)))),Gamma=Gamma_m))));

.727930005100072

1.07317073170732

5.00158254780052

N = 5 yields RL = -25.989
N = 6 yields RL = -32.9555

N := ceil(N);
# calculate new resulting Gamma_m
Gamma_m := abs(solve(sec(theta_m)=cosh((1/floor(N))*arccosh(abs(ln(ZL/ZS)/(2*Gamma)))),Gamma)[1]);

6

0.225022725352146e-1

# ONLY NEEDED IF GAMMA_M SOLUTION FAILS ABOVE
#plot([sec(theta_m),cosh((1/floor(N))*arccosh(abs(ln(ZL/ZS)/(2*Gamma))))],Gamma=0..2);

# ONLY NEEDED IF GAMMA_M SOLUTION FAILS ABOVE
#Gamma_m := abs(fsolve(sec(theta_m)=cosh((1/floor(N))*arccosh(abs(ln(ZL/ZS)/(2*Gamma)))),Gamma=0..2));

theta_m := solve(sec(theta)=cosh((1/N)*arccosh(abs(ln(ZL/ZS)/(2*Gamma_m)))),theta);

.727930005100076

G := 'G':
n := 0:
Gamma :=
  `if`(
    type(N,odd),
    eval(2*exp(-I*N*theta)*(sum('G[n]*cos((N-2*n)*theta)','n'=0..(N-1)/2-1)+G[(N-1)/2]*cos(theta))),
    eval(2*exp(-I*N*theta)*(sum('G[n]*cos((N-2*n)*theta)','n'=0..N/2-1)+G[N/2]/2))
  ):

Tsc := combine(simplify(ChebyshevT(N,sec(theta_m)*cos(theta)),'ChebyshevT')):

Gamma/exp(-I*N*theta)=Gamma_m*Tsc:  # note Gamma_m = |A|

G := 'G':
eq1 := Gamma/exp(-I*N*theta):
eq2 := Gamma_m*Tsc:
for n from 0 to trunc(N/2)-1 do
  G[n] := solve(coeff(Gamma/exp(-I*N*theta),cos((N-2*n)*theta))=coeff(Gamma_m*Tsc,cos((N-2*n)*theta)),
              G[n]);
  eq1 := eq1-coeff(Gamma/exp(-I*N*theta),cos((N-2*n)*theta))*cos((N-2*n)*theta);
  eq2 := eq2-coeff(Gamma_m*Tsc,cos((N-2*n)*theta))*cos((N-2*n)*theta);
end do:
G[trunc(N/2)] := solve(eq1=eq2,G[trunc(N/2)]):
for n from trunc(N/2)+1 to N do
  G[n] := G[N-n];
end do:
'Gamma[n]' = seq(G[n],n=0..N);

Gamma[n] = (0.649877859143785e-1, .172604292049981, .287211372499400, .336687460192378, .287211372499400, .172604292049981, 0.649877859143785e-1)

Z := 'Z':
#Z[1] := solve(ln(Z[1])=ln(ZS)+2*G[0],Z[1]):
Z[1] := evalf(exp(ln(ZS)+2*G[0])):
seq(assign(Z[n],solve(ln(Z[n])=ln(Z[n-1])+2*G[n-1],Z[n])),n=2..N);
'Z[n]' = seq(Z[n],n=1..N);
# GETS Z1 CLOSE, BUT ZN IS BAD <-- 9/15/09 - is this comment old or still valid?

Z[n] = (3.55875176303013, 5.02596984406875, 8.92664999631936, 17.5037668178350, 31.0885271594686, 43.9058440724061)

# tuning limits for optimizer
chvZtunlim := [sqrt(ZS*Z[1]),seq(sqrt(Z[n]*Z[n+1]),n=1..N-1),sqrt(Z[N]*ZL)];

[3.33483121903781, 4.22920548608316, 6.69813956931755, 12.5000000000001, 23.3273729791691, 36.9454736862906, 46.8539454434768]

Zin := ZL:
Zin_save := ZL:
for i from 0 to N-1 do
  Zin := simplify(Z[N-i]*(Zin_save+I*Z[N-i]*tan(E))/(Z[N-i]+I*Zin_save*tan(E))):
  Zin_save := Zin: # used for iteration when just plugging Zin back in doesn't work
end do:

#Zin; # should be univariate function of E

E := (Pi/2)*(f/f0): # assuming vp has no strong f-dependence
Gamma := abs((Zin-ZS)/(Zin+ZS)):

# THIS STEP COMMENTED OUT BECAUSE IT TAKES EXTREMELY LONG TO COMPUTE
# ATTEMPTING TO CALCULATE AND PLOT RESPONSE OF FILTERS WITH N > 4 MAY
# CAN SOMETIMES CAUSE MAPLE TO CRASH OR RUN OUT OF MEMORY
#if f0*(1-BW) < 0 then
#  plot(20*log10(Gamma),f=0..f0*(1+BW),y=-50..0,labels=["freq [GHz]","Gamma [dB]"],labeldirections=[horizontal,vertical])
#else
#  plot(20*log10(Gamma),f=f0*(1-BW)..f0*(1+BW),y=-50..0,labels=["freq [GHz]","Gamma[dB]"],labeldirections=[horizontal,vertical])
#end if;

 

Optimization

 

# ADDED 3/28/22
# assumes Lstub = Pi/2
with(Optimization):
ZSp := (1/ZS+1/Zstub*(1+exp(-2*j*Pi/2*f/f0))/(1-exp(-2*j*Pi/2*f/f0)))^(-1); # create parallel connection of stub and inputs

1/(8/25+(1+exp(-(.243902439024390*I)*Pi*f))/(Zstub*(1-exp(-(.243902439024390*I)*Pi*f))))

# INSERT STUB INTO CURRENT TRANSFORMER STRUCTURE
Zin := ZSp:
Zin_save := ZSp:
for i from 1 to N do
  Zin := simplify(Z[i]*(Zin_save+I*Z[i]*tan(E))/(Z[i]+I*Zin_save*tan(E))):
  Zin_save := Zin: # used for iteration when just plugging Zin back in doesn't work
#  print(evalf(subs(Zstub=Z[2],f=f1,Zin)));
end do:
#Zin;

# NEW GAMMA, WITH STUB
Gamma := abs((Zin-ZL)/(Zin+ZL)):

#plot(abs(subs(Zstub=Z[2],Zin)),f=f1..f2);

# QUICK VISUAL ON HOW IT LOOKS, SHOULD ZERO AT MID-BAND
if f0*(1-BW) < 0 then
  plot([seq(subs(Zstub=Z[k],20*log10(Gamma)),k=1..N)],f=0..f0*(1+BW),labels=["freq [GHz]","Gamma [dB]"],
    labeldirections=[horizontal,vertical],legend=[seq(k,k=1..N)])
else
  plot([seq(subs(Zstub=Z[k],20*log10(Gamma)),k=1..N)],f=f0*(1-BW)..f0*(1+BW),labels=["freq [GHz]","Gamma [dB]"],
    labeldirections=[horizontal,vertical],legend=[seq(k,k=1..N)])
end if;

# original check of internal optimization
NLPSolve(subs(Zstub=Z[2],Gamma),f=f1..f2,maximize);

[HFloat(0.37527715549854496), [f = HFloat(1.9000000131130221)]]

Gammaproc := unapply(Gamma,[f,Zstub],numeric):

bandpeak := proc(p)

  local f, sol:
  if not p::numeric then return 'procname'(p); end if;

  sol := NLPSolve(Gammaproc(f,p),f=f1..f2,
                        #'evaluationlimit'=100,
                        'method'=':-branchandbound','maximize'
                        );
  if sol[1]::numeric then

    return sol[1];
  else
    return Float(undefined);
  end if;

end proc:

bandpeak(5);

.376658504439654651

sol := NLPSolve(bandpeak(Zstub), Zstub=Z[1]..Z[N],
                      method=sqp, initialpoint = [Zstub=Z[2]]);

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

Warning, limiting number of evaluations reached

[.111740781972700001, [Zstub = HFloat(43.9058440724061)]]

bandpeak(eval(Zstub,sol[2]));

.111740781972700085

temp := NLPSolve(Gammaproc(f,eval(Zstub,sol[2])),f=f1..f2,
                         'method'=':-branchandbound','maximize',
                         'evaluationlimit'=100);

[.111740781972700085, [f = HFloat(1.9)]]

plot(bandpeak, Z[1]..Z[N], adaptive=false,
      view=0.1..0.5, gridlines=true);

optpt := eval([f,Zstub],[sol[2][],temp[2][]]);

[HFloat(1.9), HFloat(43.9058440724061)]

display(
  pointplot3d([[optpt[], Gammaproc(optpt[])]], color=red,
                   symbol=solidcircle,symbolsize=30),
  plot3d(Gammaproc(f,Zstub),f=f1..f2,Zstub=Z[1]..Z[N],
                             color=grey),
  axes=box, labels=[f,Zstub,"Gamma"]);

 

 

Download transformer_with_short_(Primes)_ac.mw

First 70 71 72 73 74 75 76 Last Page 72 of 336