tomleslie

13876 Reputation

20 Badges

15 years, 164 days

MaplePrimes Activity


These are answers submitted by tomleslie

Flip the entries in the type definition, as in

restart;
TypeTools:-AddType('type_1', `&*`(identical(x), integer));
type(3*x,type_1);

restart;
TypeTools:-AddType('type_1', '`*`'({identical(x), integer}));
type(3*x,type_1);

which returns

false
true

If I examine your PDE system, then I conclude that you need six boundary conditions. When I select the six appropriate boundary conditions from the list which you have, then the system can be solved. See the attached worksheet ( pdeProb.mw )

The only remeing problem is that if one increases the range of time over which a solution is required then convergence problems crop up. I'll leave addressing this until you confirm that the attached worksheet is OK - so far,

It states (check the highlighted part)

Error, (in fproc) unable to store '732.092516359893580*1[f]/1[nf]' when datatype=float[8]

You get this because you are using k[f] and k[nf] as parameters in your equations. You are also using 'k' as a loop index, so when the loop index is '1 k[f] becomes 1[f] and k[nf] becomes 1[nf]. But of course since 'f' and 'nf' are nowhere defined. there is no way of evaluating these to floating point values. Interestingly you have the quantities knf and kf in the list called params.

I made the wild guess that everywhere you use k[nf] and k[f], you actually mean knf and kf (from the params list). If I make these changes then the worksheet attached ( odeProb.mw ) executes without problems, although the answer from every loop iteration appears to be the same. I haven't attempted to diagnose why this happens

 

 

 

if I guess a value for alpha - I picked 0.1, then the code

  plot( [ seq( eval(c__3*t^(alpha-1)/GAMMA(alpha), alpha=0.1),
               c__3=1..9,2
             )
        ],
        t=1..10,
        legend=[ seq( typeset(c__3=j),
                      j=1..9,2
                    )
               ]
      );

produces the plot

In order to solve/plot your eq1 - you are going to have to provide a lot of further information - eg values for 'r' 'sigma' etc , etc, as well as precisely what you want to plot say u(x,t) maybe

Please do this in the form of a Maple worksheet which yu can upload using the big green up-arrow in the Mapleprimes toolbar

Somehow you have managed to enter you integrals in "inert" form, so you will not be able to evaluate these unless you "force" it with a value() command. However when you do this, all of them end up as either 0 or 1, and your final result is a tautology - but at least that means it is true!

See the attached

badSolve.mw

that there is a sensible way to grab the information you want using the PlanePlot() command. Although not shown in your "screen grab" you have obviously set infolevel=1 (or higher) whihc has the effect of showing some "intermediate" values which the PlanePlot() command has to compute in order to produce the final plot.

So I'm going to assume that all you really want is an orthonormal basis for you plane (and maybe the vector normal to the plane). The code below (and in the attached worksheet planeProb.mw) does this for specific values of your parameters alph||29: and gam||29, which are unspcified in your original post. Change these as desired

  restart:
  alias(VC=VectorCalculus, LA=LinearAlgebra, PL=plots):
#
# Set some arbitrary values for the parameters
# alph||29 and gam||29. OP does not supply these
# so just change to any desired values
#
  alph||29:=1.0:
  gam||29:=-1.0:
#
# Define the plane
#
  pl:= expand(z+tan(-Pi/2+alph||29)*(x*sin(-Pi/2+gam||29)-y*cos(-Pi/2+gam||29))=0):
#
# Need to obtain two independent vectors in the plane
#
# For any plane given by an equation of the form
#
#  A*x+B*y+C*z=0
#
# the vectors given by the components
#
#  [B/A, -1, 0] and
#  [C/A, 0, -1] are
#
#  1. guaranteed to lie in the plane, and
#  2. are not collinear
#
# Note that these are just two possibilities from an
# infinite number - so this choice is pretty arbitrary
# (and will fail if A=0 - a possibility I'm ignoring
# here). So generate these two vectors
#
  v1:= Vector( [ coeff(lhs(pl), y)/coeff(lhs(pl), x),
                 -1,
                  0
               ]
             ):
  v2:= Vector( [ coeff(lhs(pl), z)/coeff(lhs(pl), x),
                 0,
                 -1
               ]
             ):
#
# Use the Gram-Schmidt process to generate an orthonormal
# basis from the two vectors above
#
  v1n, v2n:=VC:-Vector~
                ( LA:-GramSchmidt
                      ( [v1,v2],
                        normalized=true
                      )
                )[]:
#
# Generate a (normalized) vector which is orthogonal to
# the above basis. This is a "normal" to the plane
#
  v3n:= VC:-Normalize
            ( VC:-CrossProduct
                  ( v1n,
                    v2n
                  )
            ):
#
# Plot the vectors v1n, v2n and v3n as well as the plane.
# Many options avialable for both of these plots, of which
# I have used a pretty minimal combination
#
  PL:-display
      ( VC:-PlotVector
            ( [v1n, v2n, v3n],
              color=[blue, blue, black]
            ),
        PL:-implicitplot3d
            ( pl,
              x=-1..1,
              y=-1..1,
              z=-1..1,
              style=surface,
              color=red
            )
     );
#
# Return the relevant vectors v1n, v2n, v3n
#
  v1n;v2n;v3n;

This will produce the plot below as well as the required vectors.


 

 

 

with the points
[[0.1, 0.5], [0.7, 0.7]],

neither of these should appear in the view range specified by view = [0 .. 1, 0 .. 0.1] , because both have y-values greater than 0.1.

So what exactly is your issue??

Consider the code below

  restart;
  with(plots):
  pointplot( [ [0.1, 0.5], [0.7, 0.7]],
             symbol=solidcircle,
             symbolsize=20,
             view = [0..1, 0..1]
          );
  pointplot( [ [0.1, 0.5], [0.7, 0.7]],
             symbol=solidcircle,
             symbolsize=20,
             view = [0..1, 0..0.5]
          );
  pointplot( [ [0.1, 0.5], [0.7, 0.7]],
             symbol=solidcircle,
             symbolsize=20,
             view = [0..1, 0..0.1]
          );

The view range in the first plot will allow both points to be shown

The view range in the second plot will allow one point to be shown (because the other has a y-value outside the range 0..0.5)

The view range in the third plot will show no points ( because both have a y-value outside the range 0..0.1

which you can circumvent with

`*`(1.0,`^`(10, -Digits));
                                     -10
                       1.000000000 10   

 

in the code snip you supplied

  1. the functiom 'to_c()' is defined but never used
  2. the function 'tpn() is used but never defined

If I make the rash assumption that this is a simple typo and these are intended to be the same, as in the code (change highlighted in ref)

restart:
na := 8.069439677916595*10^5:
nb := -1.777065899098942*10^3:
nc := 1.467451715287991:
nd := -5.383733471268420*10^(-4):
ne := 7.404613067985871*10^(-8):
am := 0.77317633747818500000:
bm := -0.00626025741156560000:
cm := 0.00002185947833342660:
ap := 471.909671218139000000000000:
bp := -7.368938071612570000000000:
cp := 0.041111235018593800000000:
dp := -0.000098963929768727000000:
ep := 0.88147417256725300*10^(-7):

c_y := y -> na + nb*y + nc*y^2 + nd*y^3 + ne*y^4:
to_c := ppm -> ap + bp*ppm + cp*ppm^2 + dp*ppm^3 + ep*ppm^4:
to_y := y -> to_c(c_y(y)):
ta := x -> 0.0014*int(to_y(y), y = 2000 .. x):
plot( ta(x), x=2000..2030);

this will produce the plot

I am a research scholar in mathematics
Trust me no-one her cares who/what you are

I wish to obtain the solution given in the research article
Well, you don't show a research article, just some random equations

I am uploading the pdf containing problem.
Congratulations, you manged to upload a pdf file

In the paper they mentioned that " Shooting method" is used.
welll since yiu didn't upload the paper - no-on can be sure. But it is possible tio use the "shooting method" to convert a BVP tp an IVP. which may (or may not) be useful for a BVP one of whose boundaries is "infinity". On the other hand it may, or may not be possible to solve the original BVP without going through the complexity of converting to an IVP. In addition, this conversion process is entirely numerical so would require  numerical values for all quantities in the definition of your system. SInce you provide absolutely none of these numerical values, no-one here is going to be able to help very much

The type() command and the select() command are different! You appear to think that they are identical????

type(myExpr, mytype) tests for a match of the whole expression

select( myType, myExpr) maps 'myType' across the operands of the expression.

So in your case of an  expression such as sin(x)*cos(x) and a defintion of myType as

 ''`*`'( {'specfunc'(cos),'specfunc'(sin)})';

then type(myExpr, mytype) this will match the 'whole" expression.

However since the operands are sin(x) and cos(x), neither of these operands will match the type specification because neither cos(x) nor sin(x) will match ''`*`'( {'specfunc'(cos),'specfunc'(sin)})';

you need to be looking at a structured type. A very simple-minded way to produce the two terms in your example would be

  op( select
      ( j-> type
            ( j,
              Or( `&*`
                  ( integer,
                    symbol,
                    `+`,
                    'piecewise'(t),
                    function
                  ),
                  `&*`
                  ( integer,
                    `+`,
                    'piecewise'(t)
                  )
                )
           ),
        MyExpr
      )
    );

But I would  warn you that generalising this for other cases could be "fiddly".

For example the 'integer' type is only needed in the above because both of your desired terms are products with a leading 'minus' sign - so Maple refrds the first term in the product as '-1'. If your desired terms happen to be "positive" then their will be'integer' type. Obviously you could put another couple of terms in the 'Or()' to allow for this - but I warn you than I regard using strucured types as a bit of an "experimental science".
 

It may be "neater" to define the structure types you want separately as in

`type/myType1`:=`&*`( integer, symbol, `+`, 'piecewise', function );
`type/myType2`:=`&*`( integer, `+`, 'piecewise' );
`type/myType3`:=`&*`( symbol, `+`, 'piecewise', function );
`type/myType4`:=`&*`( `+`, 'piecewise');

  op( select
      ( j-> type
            ( j,
              Or( myType1,
                  myType2,
                  myType3,
                  myType4
                )
           ),
        MyExpr
      )
    );

 

for specific values of M, N, k.

I just "guessed" some values -change the entries in the list 'params' - change these to whatever you like

See the code below and the attachment

restart;
h := k - (k - 1)*x;
DE :=   diff(p(x), x)
      = -2*M^2
        *
        ( -exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
          + M*sinh(M*h/sqrt(N))*lambda + exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)
        )
        /
        (   exp(-M*h/sqrt(N))
          *
          (   M*exp(2*M*h/sqrt(N))*h 
             - 2*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
             - M*h + 4*exp(M*h/sqrt(N))*sqrt(N)
             - sqrt(N)*exp(2*M*h/sqrt(N))
             - sqrt(N)
          )
        );
BC := p(0) = 0, p(1) = 0;
#
# Find a solution numerically using both boundary conditons
# above as well as the values in the params list for M, N,
# k, leaving lambda as unknown
# 
# The dsolve process will produce a solution for p(x) and
# the associated value of lambda. If the vlaues in params
# are changed then the value of lambda will presumably
# change, along with the value of p(x)
#
  params := [M = 0.1, N = 0.1, k = 0.1];
  sol := dsolve(eval({BC, DE}, params), numeric);
#
# Check the value of lambda
#
  sol(0.5)[3];
#
# Plot the solution for p(x). Note that both boundary
# conditions are fulfilled
#
  plots:-odeplot(sol, [x, p(x)], x = 0 .. 1);

odeProb.mw

This produces the lambda-value

lambda = 0.0908819796701764

and the plot

 

 

 

 

 

 

because (to be honest) I can't work out from your worksheet what you are trying to achieve

If it is a simple change of variables then you should probably just use the PDETools:-dchange() command as in the code below and the attachment

  restart;
#
# The PDE system in the "old" coordinates [x__o, y__o]
#
  OldSys:= [ diff(u(x__o, y__o), x__o) + diff(v(x__o, y__o), y__o) = 0,

             u(x__o, y__o)*diff(u(x__o, y__o), x__o) + v(x__o, y__o)*diff(u(x__o, y__o), y__o)
             =
               diff(u(x__o, y__o), y__o, y__o)
             +
               6*k*diff(u(x__o, y__o), y__o)^2*diff(u(x__o, y__o), y__o, y__o),

             u(x__o, y__o)*diff(T(x__o, y__o), x__o)
             +
               v(x__o, y__o)*diff(T(x__o, y__o),y__o)
             =
               diff(T(x__o, y__o), y__o, y__o)
           ]:
#
# Define the coordinate transformations
#
  tr := [ isolate(x__n = x__o/L, x__o),
          isolate(y__n = y__o*(U__infinity*L/nu)^(1/2)/L, y__o)
        ];
#
# Perform the coordinate transformations on the PDE system
#
  NewSys:= PDETools:-dchange( tr,
                              OldSys, [x__n, y__n],
                              params = [L, nu, U__infinity]
                            );

chVar.mw

 

 

I really don't understand the question!

Syrup allows you to define a network with a sequence of statements of the form

componentName nodeNumber nodeNumber componentValue

The first charcter in the componentName determine the type of the component - so R1, R2, R3 are resistors, C1, C2, C3 are capacitors etc etc.

So for the diagram you supply the equivalent Syrup descritption would be ( assuning I haven't made any typos, because I gave up typing neltists around 1988 when schematic capture packages became available and netlist generation was automated)

  ckt:="
         V1 1 0 1
         R1 1 3 50
         L1 3 2 .96e-5
         Cp 2 0 .83e-10
         L2 2 5 500e-6
         R2 5 4 0.2
         RL 4 6 1.343
         LL 6 0 0.155e-05
         .end
       ":

Having defined the circuit in this way, you just have to determine what you want to know. In the attached I have plotted the transfer function from node 1 to node 4 and the transfer function from node 1 to node 6 in decibels with a logarithmic frequency axis. It would be trivial to change the plots to something else - but since I don't know what you want I can't?!!

The code

  with(Syrup):
  with(plots):
  ckt:="
         V1 1 0 1
         R1 1 3 50
         L1 3 2 .96e-5
         Cp 2 0 .83e-10
         L2 2 5 500e-6
         R2 5 4 0.2
         RL 4 6 1.343
         LL 6 0 0.155e-05
         .end
       ":
  sol:= Solve(ckt):
  g1:= eval(v[4]/v[1], sol):
  g2:= eval(v[6]/v[1], sol):
  semilogplot( [ 20*log10(g1),
                 20*log10(g2)
               ],
               s=1..1e09,
               color=[red,blue],
               view=[1..1e09, -150..0]
             );

produces the plot

See the attached for more detail

syrup.mw

 

 

First 45 46 47 48 49 50 51 Last Page 47 of 207