tomleslie

Don't create a new matrix...

Answered: tomleslie 13876

February 21 2018

0 1

I basically agree with JoeRiel's argument about "distinct" matrices, and also that it seems like a slightly odd thing to want to do. But if you really want to achieve this, then some variation on the following "toy" example may be applicable
with(LinearAlgebra):
M:=RandomMatrix(3,3);
eq:=M+a;
eq2:=subs( op(2,eq)=A, eq);

Get the file path correct...

Answered: tomleslie 13876

February 21 2018

0 0

and once you have the file path correct, get the file name correct: and once you have the file name correct, get the worksheet names correct.

I downloaded your Excel file to my deskstop: the file path to my desktop is "C:/Users/TomLeslie/Desktop" but yours will be different. Since the file had two extensions (why?) I renamed it from Employees.xlsx.xlsx to Employees.xlsx.

So I now know that the path to the file is given by "C:/Users/TomLeslie/Desktop" (youre will be differen tdepending on where you stored the file) and the fileName is " Employees.xlsx", so when I execute

M1:=ExcelTools:-Import( "C:/Users/TomLeslie/Desktop/Employees.xlsx", "Sheet1", "A1:B101");
M2:=ExcelTools:-Import( "C:/Users/TomLeslie/Desktop/Employees.xlsx", "Sheet2", "A1:B101");
M3:=ExcelTools:-Import( "C:/Users/TomLeslie/Desktop/Employees.xlsx", "Sheet3", "A1:B101");

in Maple, I get three matrixes: one for each "sheet" in your Excel file. Whsat is difficult?

why 'pdsolve'...

Answered: tomleslie 13876

February 20 2018

1 5

You do no have a system of partial differential equations. You have a system of ordinary differential equations, for whihc the solution method would be

restart;
{sqrt(2)*sqrt(F1(t))*diff(F2(t),t) = 0};
% union {F3(t) = 0};
dsolve(%, {F1(t),F2(t),F3(t)});

Including F3(t) as a requested solution seems a little pointless, since you have already stated that F3(t)=0. Also, rather obviously, your first quation is solved with F1(t)=0, for arbitrary F2(t)

Something like...

Answered: tomleslie 13876

February 18 2018

0 0

the attached. NB I can thnk of ways to make this fail, but covering all possible cases would be difficult if not impossible (eg for polynomials of order 5 or higher). However it will work for most simple cases

>

restart;

>

#
# Defien the starting expression
#
H:=(s-z0)/((s-p0)*(s-p1));
#
# Presumably in the above 's' is some
# kind of special variable and all
# other quantities are constants/parameters.
#
# So substitute for 's' in terms of a new
# 'special' variable 'z' and pre-existing
# or new parameters/constants, then convert
# the result to normal form
#
H:=normal(subs(s=(z-1)/(T-z), H));

(1)

>

#
# Define a procedure which will convert the above
# to the form required by the OP. Note that this
# requires analytic solutions of polynomials, so
# obviously won't work if the degree is higher than
# 4. (Might also be problematic if polynomial roots
# are complex?)
#
  tidy:= proc( expr, var )
               local nsol:= [ solve
                              ( numer(expr),
                                var
                              )
                            ],
                     dsol:= [ solve
                              ( denom(expr),
                                var
                              )
                            ],
                     n1:= mul
                          ( var-nsol[i],
                            i=1..numelems(nsol)
                          ),
                     d1:= mul
                          ( var-dsol[i],
                            i=1..numelems(dsol)
                          );
               return n1/d1;
         end proc:

>	# # Run the procedure on the function, for the 'special' # variable # # tidy(H, z)

(-T+z)*(z-(T*z0+1)/(1+z0))/((z-(T*p0+1)/(1+p0))*(z-(T*p1+1)/(1+p1)))

(2)

>

Download tidyPol.mw

It depends...

Answered: tomleslie 13876

February 18 2018

1 2

so much on whether you want to exchange numeric data, a couple of commands, or whole procedures - and this is no place to give a tutorial, which would be really long, simply becuase there are so many possible alternatives.

However if you just want to exchange a couple of simple assignments, then the simplest(?) way is probably

####In Maple

restart;
with(CodeGeneration):
writeto("C:/Users/TomLeslie/Desktop/mattest.m");
Matlab([1,2,3,4], resultname="a");
Matlab([5,6,7,8], resultname="b");
close( "C:/Users/TomLeslie/Desktop/mattest.m");
writeto(terminal);

###Then in Matlab

>> run 'C:/Users/TomLeslie/Desktop/mattest.m'
>> a
>> b

Note that you will have to change the file paths in the above to something relevant for your machine

Depends what is acceptable...

Answered: tomleslie 13876

February 18 2018

0 0

Just because you can define a function (piecewise or not) does not mean that it has an'analytic' Fourier transform.

However if you are prepared to 'sample' a function, then one can (usually) compute the corresponding DFT

I have added an execution group to the attached worksheet, to show how this can be done for your case. Note that there are many options when performing/displaying a DFT - the example I attach is fairly basic

DFT.mw

Typo?...

Answered: tomleslie 13876

February 18 2018

1 0

There are 20 elements in the list 'Y', but only 19 in the list 'Z'. Either add an extra element to the Z-list, or restrict the Y-lsit as in

ScatterPlot(Y[1 .. 19], Z)

An altennative...

Answered: tomleslie 13876

February 14 2018

0 0

partly becuase I think alternatives are a good thing, and partly becuase I think the use of the Physics package is overkill in this case, when the same calculation could be performed with the (simpler?) VectorCalculus() package: although unfortunately I may have detected a bug in the latter - see the last coupel of eecution groups

>

  restart;
#
# Use the vector calculus package
#
  with(VectorCalculus):
#
# Define the two vectors
#
  r1:= t-> PositionVector([2*t^2, 16, 10*t-12], cartesian[x,y,z]):
  r2:= t-> PositionVector([20-6*t, 4*t^2, 2*t^2], cartesian[x,y,z]):
#
# Compute the time at which "endpoints" of the
# two vectors coincide. Note that this is the time
# at which the vector given by r1(t)-r2(t) has
# zero length
#
  tColl:=RealDomain:-solve(Norm(r1(t)-r2(t)),t);
#
# Generate the vectors at this time - note that since these
# are position vectors, and for "collision" to occur, both
# must result in the same position
#
  r1(tColl), r2(tColl);
#
# Since both vectors must result in the same position in
# order for a collision to occur, the angle between them
# must be zero. You don't need Maple to see this. If you
# want to use Maple then define a function which computes
# the angle between any two vectors. The Dot Product of
# two vectors a, b is given by |a|*|b|*cos(theta), so the
# angle between a and b is given by
#
# arccos( DotPoduct(a,b)/( Norm(a)*Norm(b) ) )
#
  getAngle:= (v1, v2)-> arccos( DotProduct(v1, v2)/( Norm(v1)*Norm(v2) )):
#
# Now apply this function to the position vectors
# r1, r2, at the collision time tColl. This returns
# 0 as expected
#
  getAngle( r1(tColl), r2(tColl) );

Vector(3, {(1) = 8, (2) = 16, (3) = 8}, attributes = [positionvector, coords = cartesian[x, y, z]]), Vector(3, {(1) = 8, (2) = 16, (3) = 8}, attributes = [positionvector, coords = cartesian[x, y, z]])

(1)

>

#
# Of marginally more interest, define a function which
# returns the velocity vector of a time-dependent position
# vector, then compute the two velocity vectors associated
# with the position vectors defined above, and in particular
# the two velocities at the time of collision
#
  v:= (vec, time)->D(vec)(time):
  v(r1, t),v(r2, t);
  v(r1, tColl),v(r2, tColl);

Vector(3, {(1) = 4*t, (2) = 0, (3) = 10}), Vector(3, {(1) = -6, (2) = 8*t, (3) = 4*t})

Vector[column](%id = 18446744074373287198), Vector[column](%id = 18446744074373279742)

(2)

>	# # How get the angle between the two velocities. NB I don't # think the extra eval() stage should be necessary. This # may be a Maple bug(ette) # getAngle( eval( v(r1, tColl), t=tColl), eval( v(r2, tColl), t=tColl)); evalf(%*180/Pi)

(3)

>

#
# Why is the extra eval()stage in the above necessary?
# Basically when computing the differential of a
# position vector, one *ought* to get a 'rooted' vector
# with the root point being the 'end' of the position
# vector. Maple computes the components of the rooted
# vector correctly, but does not evaluate the root point.
#
# Consider the following and note the root point expression
#
About( v(r1, tColl) );

_rtable[18446744074373274310]

(4)

>

Download vectorCalc.mw

Is cheating allowed?...

Answered: tomleslie 13876

February 12 2018

0 0

I suspect that Maple has an in-built list of "known" Mersenne primes, so for example

restart;
NumberTheory:-IsMersenne(57885161);
NumberTheory:-IsMersenne(74207281);

Both return true, more or less instantaneously. The Mersenne prime with exponent 77232917 was only discovered in December 2017, so probably hasn't made it to Maple's lookup table yet.

I suspect (but can't prove) that any index less than the highest known to Maple (ie 74207281) will return false fairly quickly. Any index greater than the highest known to Maple - well you are pretty much into the realms of testing primality of huge numbers, and that gets very, very painful, very, very quickly

2D Input quirk/feature...

Answered: tomleslie 13876

February 12 2018

0 1

Because the attached just works

Vmax = RmsMS phase to neutral voltage

=

The #1 diode is forward biased and switches on when Va>Vc (Va=Vc). The voltage VD will be equal to Va-Vb until diode #2 turns on when Vc< Vb. Then VD will equal Vc-Vc.

The switching angle for the diodes is given above. If phase control is used all commutation is delayed by the angle α.

Without phase delay the average (DC) voltage applied to the laod is the average value of the 6 commutation periods over one cycle:

Emax = RmsMS phase to neutral voltage

Eq1 := Vdo = 6*(int(Van-Vbn, theta = com1 .. com2))/(2*Pi)

Vdo = 3*3^(1/2)*2^(1/2)*Erms/Pi

(1)

Eq2 := Vdc = 6*(int(Van-Vbn, theta = com1+alpha .. com2+alpha))/(2*Pi)

Vdc = 3*2^(1/2)*Erms*3^(1/2)*cos(alpha)/Pi

(2)

Commutation

The switching between diodes can not occur instantaneously because the supply circuit is inductive. Current can not change instantaneously in and inductive circuit.
The transfer of load current from the out-going phase to the incomming phase is accomplshed by a period of symultaneous conduction of the two adjacent diodes the creation of a circulating current between these two phases that resembles short circuit between the two phases.

The circulating current produces equal voltage drop with opposid polarities in each phase. In the out-going phase the drop adds to the soruce voltage and in the incomming phase the drop subtracts from the source voltage. Only the drop in the incomming phase effects the DC output voltage.

During commutation intervil com1 to com2 phases a and c are shorted while Id is transfered from phase c to phase a. During this time voltage Va-Vb minus the voltage drop in a phase caused by the commutation. This voltage drop is 1/2 the Vac voltage.

(1/2)*2^(1/2)*Erms*sin(theta)-(1/2)*2^(1/2)*Erms*cos(theta+(1/6)*Pi)

(3)

and Va=Vc at π/6

`ΔE__Avg` = 6*(int(`Δe`, theta = (1/6)*Pi+alpha .. (1/6)*Pi+alpha+mu))/(2*Pi)

`ΔE__Avg` = 3*((1/2)*2^(1/2)*Erms*cos((1/6)*Pi+alpha)+(1/2)*2^(1/2)*Erms*sin((1/3)*Pi+alpha)-(1/2)*2^(1/2)*Erms*sin((1/3)*Pi+alpha+mu)-(1/2)*2^(1/2)*Erms*cos((1/6)*Pi+alpha+mu))/Pi

(4)

`ΔE__Avg` = (3/2)*Erms*2^(1/2)*(cos((1/6)*Pi+alpha)+sin((1/3)*Pi+alpha)-sin((1/3)*Pi+alpha+mu)-cos((1/6)*Pi+alpha+mu))/Pi

(5)

`ΔE__Avg` = (3/2)*2^(1/2)*Erms*3^(1/2)*cos(alpha)/Pi-(3/2)*Erms*2^(1/2)*3^(1/2)*cos(alpha)*cos(mu)/Pi+(3/2)*Erms*2^(1/2)*3^(1/2)*sin(alpha)*sin(mu)/Pi

(6)

`ΔE__Avg` = (3/2)*3^(1/2)*Erms*((-cos(mu)+1)*cos(alpha)+sin(alpha)*sin(mu))*2^(1/2)/Pi

(7)

`ΔE__Avg` = (1/2)*(-3*3^(1/2)*2^(1/2)*Erms*cos(alpha+mu)+3*2^(1/2)*Erms*3^(1/2)*cos(alpha))/Pi

(8)

`ΔE__Avg` = -(3/2)*3^(1/2)*2^(1/2)*Erms*(cos(alpha+mu)-cos(alpha))/Pi

(9)

DC output voltage

(10)

V__D = 3*2^(1/2)*Erms*3^(1/2)*cos(alpha)/Pi+(3/2)*3^(1/2)*2^(1/2)*Erms*(cos(alpha+mu)-cos(alpha))/Pi

(11)

V__D = (3/2)*3^(1/2)*2^(1/2)*Erms*(cos(alpha)+cos(alpha+mu))/Pi

(12)

Equating the intigral of the voltage drop between the phase A and C to the intigral of the voltage rise across the two phase inductances (thes intigrals represent eqivalent flux).

int(Van-Vcn, theta = (1/6)*Pi+alpha .. (1/6)*Pi+alpha+mu) = int(2*L*omega, theta = 0 .. Id)

2^(1/2)*Erms*cos((1/6)*Pi+alpha)+2^(1/2)*Erms*sin((1/3)*Pi+alpha)-2^(1/2)*Erms*cos((1/6)*Pi+alpha+mu)-2^(1/2)*Erms*sin((1/3)*Pi+alpha+mu) = 2*L*omega*Id

(13)

(14)

(15)

(16)

(17)

(18)

V__D = (3/2)*3^(1/2)*2^(1/2)*Erms*(cos(alpha)-(1/6)*(2*L*omega*Id-2^(1/2)*Erms*3^(1/2)*cos(alpha))*3^(1/2)*2^(1/2)/Erms)/Pi

(19)

V__D = 3*2^(1/2)*Erms*3^(1/2)*cos(alpha)/Pi-3*L*omega*Id/Pi

(20)

Line Currents

A phase current intervals

Commutation from C to A (from π/6+ α to π/6 + α+ μ)

Initially A phase current is zero and C phase current = Idc. Commutation begins when Va>Vc plus any delay angle. The commutation current (iac) is
equivalent to a short circuit between C and A phase through 2 x the phase reactance.

since then

iac = (1/2)*(-2^(1/2)*Erms*cos(theta)-2^(1/2)*Erms*sin(theta+(1/6)*Pi))/(L*omega)

(1.1)

At θ= π/6+ α Ia = 0 then Ia = iac(θ) - iac(θ= π/6+ α)

Ia = (1/2)*(-2^(1/2)*Erms*cos(theta)-2^(1/2)*Erms*sin(theta+(1/6)*Pi))/(L*omega)-(1/2)*(-2^(1/2)*Erms*cos((1/6)*Pi+alpha)-2^(1/2)*Erms*sin((1/3)*Pi+alpha))/(L*omega)

(1.2)

Ia = -(3/4)*2^(1/2)*Erms*cos(theta)/(L*omega)-(1/4)*2^(1/2)*Erms*sin(theta)*3^(1/2)/(L*omega)+(1/2)*2^(1/2)*Erms*3^(1/2)*cos(alpha)/(L*omega)

(1.3)

Ia = -(1/4)*2^(1/2)*((sin(theta)-2*cos(alpha))*3^(1/2)+3*cos(theta))*Erms/(L*omega)

(1.4)

(1.5)

(1.6)

(1.7)

(1.8)

(1.9)

(1.10)

Ia = -(1/4)*2^(1/2)*((sin((1/6)*Pi+alpha)-2*cos(alpha))*3^(1/2)+3*cos((1/6)*Pi+alpha))*Erms/(L*omega)

(1.11)

(1.12)

For RMS calc area under Ia^2

A1 = int(rhs(Ia = -(1/4)*2^(1/2)*((sin(theta)-2*cos(alpha))*3^(1/2)+3*cos(theta))*Erms/(L*omega))^2, theta = (1/6)*Pi+alpha .. (1/6)*Pi+alpha+mu)

(1.13)

A1 = -(3/16)*Erms^2*a1/(L^2*omega^2)

(1.14)

(1.15)

a1 = (-8*cos(mu)*sin(mu)-8*mu+16*sin(mu))*cos(alpha)^2-8*sin(alpha)*(cos(mu)-1)^2*cos(alpha)+4*cos(mu)*sin(mu)-4*mu

(1.16)

A1 = -(3/16)*Erms^2*((-8*cos(mu)*sin(mu)-8*mu+16*sin(mu))*cos(alpha)^2-8*sin(alpha)*(cos(mu)-1)^2*cos(alpha)+4*cos(mu)*sin(mu)-4*mu)/(L^2*omega^2)

(1.17)

A phase carries Idc from π/6+ α + μ to 5π/6 +α

A2 = int(Idc^2, theta = (1/6)*Pi+alpha+mu .. 5*Pi*(1/6)+alpha)

(1.18)

Commutation from A to B (from π/5+ α to 5π/6 + α+ μ)

Initially A phase current equals Idc and B phase current = 0. Commutation begins when Va<Vc plus any delay angle. The commutation current (iac) is
equivalent to a short circuit between C and A phase through 2 x the phase reactance.

since then

iba = (1/2)*(2^(1/2)*Erms*cos(theta)+2^(1/2)*Erms*cos(theta+(1/3)*Pi))/(L*omega)

(1.19)

At θ= 5π/6+ α Ia = Idc then Ia = iac(θ) - iac(θ= 5π/6+ α)+ Idc

Ia = -(1/2)*(2^(1/2)*Erms*cos(theta)+2^(1/2)*Erms*cos(theta+(1/3)*Pi))/(L*omega)+(1/2)*(2^(1/2)*Erms*cos((5/6)*Pi+alpha)+2^(1/2)*Erms*cos((7/6)*Pi+alpha))/(L*omega)+Idc

(1.20)

Ia = (-(1/2)*(cos(theta)*2^(1/2)+2^(1/2)*cos(theta+(1/3)*Pi))/(L*omega)+(1/2)*(-cos((1/6)*Pi+alpha)*2^(1/2)-sin((1/3)*Pi+alpha)*2^(1/2))/(L*omega))*Erms+Idc

(1.21)

Ia = (-(1/2)*(cos((5/6)*Pi+alpha)*2^(1/2)+2^(1/2)*cos((7/6)*Pi+alpha))/(L*omega)+(1/2)*(-cos((1/6)*Pi+alpha)*2^(1/2)-sin((1/3)*Pi+alpha)*2^(1/2))/(L*omega))*Erms+Idc

(1.22)

(1.23)

Ia = (-(1/2)*(cos((5/6)*Pi+alpha+mu)*2^(1/2)+2^(1/2)*cos((7/6)*Pi+alpha+mu))/(L*omega)+(1/2)*(-cos((1/6)*Pi+alpha)*2^(1/2)-sin((1/3)*Pi+alpha)*2^(1/2))/(L*omega))*Erms+Idc

(1.24)

(1.25)

(1.26)

For RMS calc area under Ia^2

A3 = int(rhs(Ia = (-(1/2)*(cos(theta)*2^(1/2)+2^(1/2)*cos(theta+(1/3)*Pi))/(L*omega)+(1/2)*(-cos((1/6)*Pi+alpha)*2^(1/2)-sin((1/3)*Pi+alpha)*2^(1/2))/(L*omega))*Erms+Idc)^2, theta = 5*Pi*(1/6)+alpha .. 5*Pi*(1/6)+alpha+mu)

(1.27)

A3 = (1/4)*(4*3^(1/2)*omega*Idc*Erms*((sin(mu)-mu)*cos(alpha)+sin(alpha)*(cos(mu)-1))*L*2^(1/2)+6*Erms^2*(cos(mu)*sin(mu)-2*sin(mu)+mu)*cos(alpha)^2+6*sin(alpha)*Erms^2*(cos(mu)-1)^2*cos(alpha)+4*Idc^2*L^2*mu*omega^2-3*Erms^2*cos(mu)*sin(mu)+3*Erms^2*mu)/(L^2*omega^2)

(1.28)

A3 = (1/4)*(4*3^(1/2)*omega*Idc*Erms*((sin(mu)-mu)*cos(alpha)+sin(alpha)*(cos(mu)-1))*L*2^(1/2)+6*Erms^2*(cos(mu)*sin(mu)-2*sin(mu)+mu)*cos(alpha)^2+6*sin(alpha)*Erms^2*(cos(mu)-1)^2*cos(alpha)+4*Idc^2*L^2*mu*omega^2-3*Erms^2*cos(mu)*sin(mu)+3*Erms^2*mu)/(L^2*omega^2)

(1.29)

A3 = (1/4)*(6*(cos(mu)*sin(mu)-2*sin(mu)+mu)*cos(alpha)^2+6*sin(alpha)*(cos(mu)-1)^2*cos(alpha)-3*cos(mu)*sin(mu)+3*mu)*Erms^2/(L^2*omega^2)+3^(1/2)*Idc*((sin(mu)-mu)*cos(alpha)+sin(alpha)*(cos(mu)-1))*2^(1/2)*Erms/(omega*L)+Idc^2*mu

(1.30)

RMS Current Line

(2.1)

I__RMS = sqrt(A/Pi)

I__RMS = (A/Pi)^(1/2)

(2.2)

-(3/16)*Erms^2*((-8*cos(mu)*sin(mu)-8*mu+16*sin(mu))*cos(alpha)^2-8*sin(alpha)*(cos(mu)-1)^2*cos(alpha)+4*cos(mu)*sin(mu)-4*mu)/(L^2*omega^2)+Idc^2*((2/3)*Pi-mu)+(1/4)*(6*(cos(mu)*sin(mu)-2*sin(mu)+mu)*cos(alpha)^2+6*sin(alpha)*(cos(mu)-1)^2*cos(alpha)-3*cos(mu)*sin(mu)+3*mu)*Erms^2/(L^2*omega^2)+3^(1/2)*Idc*((sin(mu)-mu)*cos(alpha)+sin(alpha)*(cos(mu)-1))*2^(1/2)*Erms/(omega*L)+Idc^2*mu

(2.3)

subs([(-8*cos(mu)*sin(mu)-8*mu+16*sin(mu))*cos(alpha)^2-8*sin(alpha)*(cos(mu)-1)^2*cos(alpha)+4*cos(mu)*sin(mu)-4*mu = a, (6*(cos(mu)*sin(mu)-2*sin(mu)+mu))*cos(alpha)^2+6*sin(alpha)*(cos(mu)-1)^2*cos(alpha)-3*cos(mu)*sin(mu)+3*mu = b, (sin(mu)-mu)*cos(alpha)+sin(alpha)*(cos(mu)-1) = c], A);

-(3/16)*Erms^2*a/(L^2*omega^2)+Idc^2*((2/3)*Pi-mu)+(1/4)*b*Erms^2/(L^2*omega^2)+3^(1/2)*Idc*c*2^(1/2)*Erms/(omega*L)+Idc^2*mu

(2.4)

Download subsProb.mw

First things first...

Answered: tomleslie 13876

February 12 2018

0 1

restart;
m := 40:
L := 5:
E :=70*(1-(y+0.5)^2)+380*(y+0.5)^2:
beta := Pi^2/L^2:
phi := add(a[n]*y^n, n = 0 .. m):
Eq := diff(phi, y$2)+(diff(E, y))*(diff(phi, y))/E+((diff(E, y$2))/E-((diff(E, y))/E)^2)*phi-2*beta*(1.3)*(phi-1):
st := [seq(coeftayl(Eq, y = 0, j), j = 0 .. m-2)]:
numEqs:=numelems(st);
numVars:=numelems(indets(st));

indicates that you have 39 equations in 41 unknowns - is this deliberate?

If it is deliberate then any subsequent solve() process will probably require to know which of the 39 variables you wish to solve for (and by inference, which of the two variables are 'free')

Convert to a "shooting" problem...

Answered: tomleslie 13876

February 12 2018

0 0

for which the "quick and dirty" solution would be something like

>

restart;
f:= alpha -> rhs( dsolve( [ diff(y(x), x$3)+diff(y(x),x$2)+y(x)=1,
                            y(0)=0,
                            D(y)(0)=alpha,
                            y(2)=1
                          ],
                          numeric
                        )(1.0)[2]);
dsol:= dsolve( [ diff(y(x), x$3)+diff(y(x), x$2)+y(x)=1,
                 y(0)=0,
                 D(y)(0)=fsolve(f),
                 y(2)=1
               ],
               numeric
             );
plots:-odeplot(dsol, 0..2);

proc (alpha) options operator, arrow; rhs((dsolve([diff(y(x), `$`(x, 3))+diff(y(x), `$`(x, 2))+y(x) = 1, y(0) = 0, (D(y))(0) = alpha, y(2) = 1], numeric))(1.0)[2]) end proc

$proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .2826168901444572, (3) = .5873070283757869, (4) = .8880663580130357, (5) = 1.1808504094981938, (6) = 1.4629933438259606, (7) = 1.7347212350044052, (8) = 2.0}, datatype = float[8], order = C_order); Y := Matrix(8, 3, {(1, 1) = .0, (1, 2) = -.4955860224, (1, 3) = .9603222325830569, (2, 1) = -.10145930955515191, (2, 2) = -.22123989371797662, (2, 3) = .9847554156360278, (3, 1) = -.12259359308266558, (3, 2) = 0.8434439320665829e-1, (3, 3) = 1.0203588945637323, (4, 1) = -0.50693766919627146e-1, (4, 2) = .39473497671368013, (4, 3) = 1.0391272932606694, (5, 1) = .10932981097212124, (5, 2) = .6975882755164673, (5, 3) = 1.022639104637302, (6, 1) = .3462030819415321, (6, 2) = .9786368092012532, (6, 3) = .9613367509425015, (7, 1) = .6464325836130244, (7, 2) = 1.2262024713463255, (7, 3) = .8521602972370497, (8, 1) = 1.0, (8, 2) = 1.4324360873364448, (8, 3) = .6940344221700293}, datatype = float[8], order = C_order); YP := Matrix(8, 3, {(1, 1) = -.4955860224, (1, 2) = .9603222325830569, (1, 3) = 0.39677767416943066e-1, (2, 1) = -.22123989371797662, (2, 2) = .9847554156360278, (2, 3) = .11670389391912417, (3, 1) = 0.8434439320665829e-1, (3, 2) = 1.0203588945637323, (3, 3) = .1022346985189333, (4, 1) = .39473497671368013, (4, 2) = 1.0391272932606694, (4, 3) = 0.11566473658957688e-1, (5, 1) = .6975882755164673, (5, 2) = 1.022639104637302, (5, 3) = -.13196891560942325, (6, 1) = .9786368092012532, (6, 2) = .9613367509425015, (6, 3) = -.3075398328840335, (7, 1) = 1.2262024713463255, (7, 2) = .8521602972370497, (7, 3) = -.49859288085007414, (8, 1) = 1.4324360873364448, (8, 2) = .6940344221700293, (8, 3) = -.6940344221700294}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .2826168901444572, (3) = .5873070283757869, (4) = .8880663580130357, (5) = 1.1808504094981938, (6) = 1.4629933438259606, (7) = 1.7347212350044052, (8) = 2.0}, datatype = float[8], order = C_order); Y := Matrix(8, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.953480294877803e-8, (2, 1) = 0.15600067486957252e-7, (2, 2) = -0.26102275883375104e-7, (2, 3) = 0.4680234128992737e-7, (3, 1) = 0.1624401045921512e-7, (3, 2) = -0.35066768178446943e-7, (3, 3) = 0.6206342268545329e-7, (4, 1) = 0.5473347330689072e-8, (4, 2) = -0.24449526320775286e-7, (4, 3) = 0.5482230913607099e-7, (5, 1) = -0.4712304738762939e-8, (5, 2) = -0.7958716847032882e-8, (5, 3) = 0.4412705245108721e-7, (6, 1) = -0.9377198271661575e-8, (6, 2) = 0.8772755330224025e-8, (6, 3) = 0.35471727418284805e-7, (7, 1) = -0.7625987607623473e-8, (7, 2) = 0.23920871135243922e-7, (7, 3) = 0.2872636032436993e-7, (8, 1) = .0, (8, 2) = 0.3730857081916653e-7, (8, 3) = 0.2233625805051706e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(6.206342268545329e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 8, [y(x), diff(y(x), x), diff(diff(y(x), x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 3, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[y(x), diff(y(x), x), diff(diff(y(x), x), x)]'[i] = yout[i], i = 1 .. 3)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(6.206342268545329e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 3, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(8, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [x, y(x), diff(y(x), x), diff(diff(y(x), x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[y(x), diff(y(x), x), diff(diff(y(x), x), x)]'[i] = res[i+1], i = 1 .. 3)] catch: error end try end proc$

>

Download shoot.mw

Made some further changes...

Answered: tomleslie 13876

February 12 2018

0 0

to Kitonum's worksheet

These are based on "semi-educated" guesses as to what the OP originally intended. Obviously these guesses may be wrong, but at least the resulting worksheet now competes with no errors

odeProb.mw

Use indets()...

Answered: tomleslie 13876

February 09 2018

0 0

as in

eq:=2*x+3*y=0;
indets(eq);

There are various options which can be supplied to the indets() command, depending exactly on what you are looking for

If I understand correctly...

Answered: tomleslie 13876

February 09 2018

0 1

and assuming I haven't made a typo somewhere, then the attached *seems* to work

fplot.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

Don't create a new matrix...

Get the file path correct...

why 'pdsolve'...

Something like...

It depends...

Depends what is acceptable...

Typo?...

An altennative...

Is cheating allowed?...

2D Input quirk/feature...

First things first...

Convert to a "shooting" problem...

Made some further changes...

Use indets()...

If I understand correctly...

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