## 5363 Reputation

9 years, 302 days

## Not an expert with Physics[]...

There are some notational niceties associated with using the Physics[Vectors] package.

In particullar, because both the Divergence() and Laplacian() commands expect their arguments to be vectors, you have to indicate that f(x,y,z) is a vector-valued function. In your worksheet, f(x,y,z) is a "simple" scalar function.

The simplest(?) way to denote a vector function within the Physics[Vectors] package is to use a postfix underscore (ie _) on the function name, so f_(x,y,z), denotes a vector function.

Thus

restart;
with(Physics[Vectors]):
Setup(mathematicalnotation=true):

ought to get you started. Whether or not your equation can be solved will ( I assume) depend on the nature of f_()

## Close to Rouben's but not quite the same...

For the solution to this problem, I came up with

restart;
sys := { diff(x(t), t) = y(t),
diff(y(t), t) = -(1+(1/100)*t)^2*x(t),
x(0) = -1,
y(0) = 2}:
dsn := dsolve( sys, numeric, output=listprocedure):
E:=  (a, b, w)->0.5*(a^2+w^2*b^2):
II:= (a, b, w)-> E(a, b, w)/w:
plots[odeplot](dsn, [  [t, E(  diff(x(t),t), x(t), 1+t/100), color=red],
[t, II( diff(x(t),t), x(t), 1+t/100), color=blue]
],
t=0..10
);

The resulting graphs disagree slighlty with Rouben's. I think that the discrepancy is due to the fact that Rouben's solution defines the function E() as

E := 1/2*(y(t)^2 + w^2*x(t)^2);

whereas the OP defined it as (my Highlight)

E(t)=0.5*(diff(x(t),t)^2+w^2*x(t)^2)

I have followed the OP's definition

## Why?...

Like Acer has already said, I have no idea why you would want to do this in the way which you specify.

You claim that

I don't want to set the result to mat in the loop to avoid memory issues while the terms in mat are calculated.

But if you have a "memory issue" then this will also show up on the read.

If you really, really want to do this (and I simply do not understand why??), then something like the following ought to work

#
# Toy example for write: Note that this does *not*
# assume that what is being written is a matrix: it
# simply writes each entry on a new line
#
# The target file (fname) is a location specific to
# my machine - OP will have to use something which works
# on his/her machine
#
restart:
fname:="J:/Users/TomLeslie/ThisIsATest":
for i from 1 by 1 to 2 do
for j from 1 by 1 to 3 do
fprintf( fname, "%a\n", x^i+2*y^j):
od;
od:
close(fname);
#
# Toy example for read. This just does a sequential
# read of the entries in the file and allocates
# them to a matrix entry, essentially based on the
# position in the specified file
#
restart;
fname:="J:/Users/TomLeslie/ThisIsATest":
A:=Matrix(2,3):
for i from 1 by 1 to 2 do
for j from 1 by 1 to 3 do
A[i,j]:=fscanf( fname, "%a")[]:
od;
od:
A;

## Many ways to do thiis...

One way is to use the objects from the geom3d package, as in

restart;
with(geom3d):
#
# Define required objects, using a small
# sphere at the intersection point
#
plane(p1,x-2*y+z=0, [x, y, z]):
plane(p2,2*y-8*z=8, [x, y, z]):
plane(p3, -4*x+5*y+9*z=-9, [x, y, z]):
sphere(s, [intersection(pt, p1, p2, p3), 1]):
#
# return cordinates of intersection point
#
ip:= coordinates(pt);
#
# draw planes and point
#
draw( [ p1(color=red),
p2(color=blue),
p3(color=green),
s(color=black)
],
view=[ seq( j-10..j+10,
j in ip
)
],
axes=boxed
);

## No exact solution?!...

Gave your problem a quick try and couldn't generate an "exact", ie symbolic,  solution.

Solving the ODE numerically is OK over the range 0..1, but solver hits a singularity at ~1.008

See attached

deProb.mw

What makes you thin an "exact" solution exists?

## Well -quite a few things...

DEplot is a command within the DETools package which must be loaded, so you either need

with(DETools);
DEplot(sys, [a(t), b(t), c(t)], t = 0 .. 2, a = -15 .. 15, b = -15 .. 15, c = -15 .. 15, color = magnitude, title = `Stable Limit Cycles`, arrows = curve, dirfield = 800, axes = none);

or

DETools[DEplot](sys, [a(t), b(t), c(t)], t = 0 .. 2, a = -15 .. 15, b = -15 .. 15, c = -15 .. 15, color = magnitude, title = `Stable Limit Cycles`, arrows = curve, dirfield = 800, axes = none);

But, neither of these will work - you will get the error message

Error, (in DEtools/DEplot) direction fields may only be produced for systems with two or fewer dependent variables.

and you have three dependent variables, a(t), b(t), and c(t). So this is never going to work.

Similarly odeplot() is a command within the plots() package, so you need either

with(plots);
odeplot(sys, [t, a(t), b(t), c(t)], -4 .. 4, color = orange);

or

plots[odeplot](sys, [t, a(t), b(t), c(t)], -4 .. 4, color = orange);

But, neither of these will work, because you have not specified in your dsolve() command that a numeric solution is being sought - so you get the error

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

Now if you want to generate analytic solutions, for plotting with subsequent commands, all you need is

dsolve(sys);

which will return

[{a(t) = a(t)}, {b(t) = b(t)}, {c(t) = _C1}]

In other words your differential equations are satisfied whenever c(t) is an arbitrary constant, and a(t), b(t) are arbitrary functions. I suspect this is not the answer you expect. Are you sure that the original three ODEs are correct??

## It does work!...

This will compute all the i+n values - but they are not being stored anywhere or assigned to anything so they "disappear" as fast as they are created. Compare with

restart;
ans:=Array(0..10,0..10):
for i from 0 to 10 do
for n from 0 to 10 do
ans[i,n]:= i+n
end do;
end do:
ans;

## hmmm - possible bug...

According to the Maple help, allowable options for the plottools[parallelepiped] command should be all of those applicable in plot3D/options - of which shading is definitely one. Other options such as color and transparency are accepted. Although one gets a warning messaage about 'shading' not being used by plottools, it is obvious from the reulting plot that some shading option has in fact been applied.

Contemplate the following

with(plottools):
with(plots):
p1:=parallelepiped([0, 0, 1], [0, 2, 2], [2, 2, 2], color=red, transparency=0.5):
p2:=parallelepiped([0, 0, 1], [0, 1, 1], [1, 1, 1], color=blue, transparency=0.5):
display([p1,p2]);
p1:=parallelepiped([0, 0, 1], [0, 2, 2], [2, 2, 2], shading=xyx, transparency=0.5):
p2:=parallelepiped([0, 0, 1], [0, 1, 1], [1, 1, 1], shading=xyz, transparency=0.5):
display([p1,p2]);

Could be a bug????

## Quick and dirty...

If you are happy with quick and dirty then

restart;
f:=x->(7-x)*sin(x^2-7):
f1:=D(f):
llim:=-2.4:
ulim:=2.4:
sol:=[RootFinding[NextZero](f1, llim)];
while true do
rt:= RootFinding[NextZero](f1, sol[-1]):
if      rt < ulim
then sol:=[sol[], rt];
else break;
fi;
od:
sol;

ought to work. You can change the start and end values of the search by fiddling with llim and ulim.

(The above is not very efficient because lists in Maple are not mutable, so the above involves a lot of unnecessary list manipulations, but if the number of entries is small, then you are probably not going to notice any time penalty)

## Just following your instructions...

restart;
#
# Starting with the equation
#
H := 1/(2*m)*diff(f(q,P,t), q)^2 = -diff(f(q,P,t),t);
#
# Make the substitution
#
#   f(q,P,t) = f1(q) + f2(t)
#
# in H and reassign the output to H
#
#
H:=subs( f(q,P,t) = f1(q) + f2(t), H);
#
# Set both sides of this equation to the same
# constant and treat as a set of two differntial
# equations - and solve this ODE system
#
sol:=dsolve({rhs(H)=E,lhs(H)=E});
#
# Select the positive solution and compute
# the combination f1+f2
#
S:=rhs(sol[1][1]+sol[1][2]);
#
# perform the requested differentiations
#
p:=diff(S,q);
Q:=diff(S,E);

## Hmmmm...

Like the manual says.

The ShowSolution command is used to show the solution steps for a Calculus1 problem, that is, a limit, differentiation or integration problem such as can be expected to be encountered in a single-variable calculus course.

so it may not be appropriate for "polynomial factoring questions and linear equation questions":-(

The second problem you might have using it, is that you have to ensure that you are preventing Maple from doing its automatic simplification: so for example

Student[Calculus1][ShowSolution]( diff(x^2*sin(x),x));

will produce an error, because the argument will be simplified (ie the diffentiation will be carried out) before being presented to ShowSolution(). This can be avoided by using uneval quotes, which postpones the simplification, so

Student[Calculus1][ShowSolution]( 'diff(x^2*sin(x),x)');

will produce a step-by-step output

For "polynomial factoring questions and linear equation questions" it may be better to take a look at the commands

Student[Basics][ExpandSteps] (), and

Student[Basics][LinearSolveSteps]

The first of these show the step-by-step guide to "simplifying" expressions involving polynomials. As an example, consider the difference between the outputs to the following two commands

simplify(expand((a^2-1)/(a/3+1/3)));
Student[Basics][ExpandSteps]( (a^2-1)/(a/3+1/3) );

The second command shows the steps involved in "solving" an equation for a variable (provided that the equation s linear in the desired variable). As an example of the second, consider the difference between the outputs to the following two commands

solve( (x+1)/(2*y*z) = 4*y^2/z + 3*x/y, x );
Student[Basics][LinearSolveSteps]( (x+1)/(2*y*z) = 4*y^2/z + 3*x/y, x );

## Lots of variables...

Depending on several factors such as,

the Maple interface you are using (1-D math, 2-D math, whatever)

the results of attempting a cut/paste operation can be pretty variable.

One thing which pretty much always works is to

create a Maple worksheet

upload that worksheet, using the big green up-arrow icon in the MaplePrimes entry window

## Errrr...

Most of this question seems to have gone missing, but something like the following will perform the basic function of plotting a function, and then plotting secants from a given anchot point

with(plots):
with(plottools):
f:=x->x^3-x^2-4*x+10: # define the function
p1:=[0, f(0)]: # the "anchor" point
p2:=[4, f(4)]: # furthest point for which the secant is requred
#
# generate a sequence of plots comprising the graph of the function
# combined with a secant

l1:= [ seq
( display( [ f1, line( p1, [t,f(t)],
color=red
)
]
),
t = p2[1]..p1[1], -1/10
)
]:
#
# display the sequence of secants
#
display(l1, insequence=true);

## Essentially an optimisation problem...

Your original code finds a numeric soltuion for the differential equation, for any supplied value of the parameter A. However you wish to specify that with theta=5/1000, determine a value for A, which makes r(theta) as close as possible to 10000 - at least that is how I understand the question.

Essentially this is an optimisation problem: if you think of r(theta), as also depending on the value of A, so very loosely, one could write it as r(theta, A), then you are looking for the situation where (with theta=5/1000)

( r(5/1000, A) -10000 )^2

is minimised - (you want it to be zero)

The issue then becomes one of defining a suitable target function for the optimisation and using NLPSolve to provide the answer (Probably could have used DirectSearch - but I didn't, cos NLPSolve worked)

See the attached code, (excuse the 1-D math input - I'm a bit old-fashioned)

dsolPar.mw

which takes about 5sec on my machine to produce the answer.

199.961306129991433161445051494

More (or fewer) digits could be obtained by fiddling with the Digits setting in the worksheet

## Just use subs()...

eg

restart;
ode:=diff(x(t),t)=x(t)*sin(t);
subs(_C1=A, dsolve(ode));

 First 107 108 109 110 111 112 113 Last Page 109 of 119
﻿