mmcdara

1247 Reputation

13 Badges

4 years, 72 days

MaplePrimes Activity


These are questions asked by mmcdara

Hi, 

Could anyone help me to numerically solve this ode?
I've tried almost all the methods Maple proposes, trying to adjust stepsizes, tolerances and so on;, always without success.

I give also the exact solution of this ode in order to compare the numerical solution to.

Thanks in advance
 

restart

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

with(plots):

# Source term

F := t*(-600*t/(100*t^2+1)^2+80000*t^3/(100*t^2+1)^3)/(100*t^2+1)-(1/(100*t^2+1)-200*t^2/(100*t^2+1)^2)/(1+(1/(100*t^2+1)-200*t^2/(100*t^2+1)^2)^2);

plot(F, t=0..0.5);

t*(-600*t/(100*t^2+1)^2+80000*t^3/(100*t^2+1)^3)/(100*t^2+1)-(1/(100*t^2+1)-200*t^2/(100*t^2+1)^2)/(1+(1/(100*t^2+1)-200*t^2/(100*t^2+1)^2)^2)

 

 

# Ode

ode := X(t)*diff(X(t), t$2)-diff(X(t),t)/(1+diff(X(t),t)^2) - 'F'

X(t)*(diff(diff(X(t), t), t))-(diff(X(t), t))/(1+(diff(X(t), t))^2)-F

(2)

# Initial conditions

ics := X(0) = 0, D(X)(0) = 1

X(0) = 0, (D(X))(0) = 1

(3)

# I used alot methods with allways either failure or either a HFloat(undefined)

printf("rkf45\n");
sol := dsolve({ode, ics}, numeric):
sol(1e-8);

printf("\n\nrosenbrock\n");
sol := dsolve({ode, ics}, numeric, method=rosenbrock):
sol(1e-8);

printf("\n\ngear\n");
sol := dsolve({ode, ics}, numeric, method=gear):
sol(1e-8);

printf("\n\ngear\n");
sol := dsolve({ode, ics}, numeric, method=classical[heunform]):
sol(1e-8);

rkf45

Error, (in sol) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

 



rosenbrock

Error, (in dsolve/numeric/SC/firststep) unable to evaluate the partial derivatives of f(x,y) for stiff solution

 

Error, (in sol) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

 



gear

 

[t = 0.1e-7, X(t) = HFloat(HFloat(undefined)), diff(X(t), t) = HFloat(HFloat(undefined))]

 



gear

 

[t = 0.1e-7, X(t) = HFloat(HFloat(undefined)), diff(X(t), t) = HFloat(HFloat(undefined))]

(4)

# The solution must be this one

U := t -> t/((t*10)^2+1)

proc (t) options operator, arrow; t/(100*t^2+1) end proc

(5)

# Check ode and ics

eval(ode, X(t)=U(t));
U(0);
D(U)(0);

0

 

0

 

1

(6)

# Plots when a solution is obtained

display(
  plot(U(t), t=0..1, color=blue),
  odeplot(sol, [t, x(t)], t=0..1, color=red, linestyle=3)
);


 

Download Unsuccessful_dsolve.mw

 

Hi, 

 

Let S(N) the set of all N by N matrices defined this way:

  • each element of a matrix M in S(N) is an integer number between 1 and included N^2
  • all the elements of M are different

For instance the matrix M = < <1, 2>|<3, 4> > belongs to S(2).

I'm interesting in finding the number of singular matrices of S(N), and more reasonnably of S(N <=3).
It's easy to verify that no matrix in S(2) is singular.
For S(3) now: as S(3) contains only 9! = 362880 elements a brute force approach can still be used. It showed that 2736 matrices were singular (about 0.75%).
But I wonderded if a more elegant approach could be used?
In the attached file I wrote all the relations elements of S(2) (and next S(3)) must verify and solved these equations for integer solutions (I only accounted for singular matrices . The case S(2) is tractable, but not S(3) (at least on my computer).

So my question: do you have any idea of some method to tackle this problem, or are you aware of any theoritical results about this issue?

PS: of course a statistical approach in which elements of S(N) would be generated using random permutations of [$1..N^2] is still possible to get a crude approximation of the number of singular matrices, but I'm not interested in this kind of approach.

Thanks in advance


 

restart:

alias(det = LinearAlgebra[Determinant])

det

(1)

 

Brute Force

 

S    := 0:
p    := 3:
PERM := combinat:-permute(p^2):
for perm in PERM do
  M := Matrix(p, p, perm):
  if det(M)=0 then S := S+1; end if:
end do:
S;
evalf(S/9!);

2736

 

0.7539682540e-2

(2)

 

A non systematic approach

Case of S(2)

 

M := Matrix(2, 2, symbol=m):
iM := {indices(M)}:

# set of all relations that define the elements of S(2)

rels :=
  {
    det(M) = 0,

    # each term is equal to an integer between 1 and 4 included

    mul((M[1, 1]-k), k=1..4)=0,
    mul((M[1, 2]-k), k=1..4)=0,
    mul((M[2, 1]-k), k=1..4)=0,
    mul((M[2, 2]-k), k=1..4)=0,

    # the sum of all the terms is equal to 10

    # M[1, 1]+M[1, 2]+M[2, 1]+M[2, 2]=10,

    # all the terms are different

    seq( mul(seq((M[op(im)]-M[op(ij)]), ij in iM minus {im})) <> 0, im in iM)
  }

{(m[1, 1]-1)*(m[1, 1]-2)*(m[1, 1]-3)*(m[1, 1]-4) = 0, (m[1, 2]-1)*(m[1, 2]-2)*(m[1, 2]-3)*(m[1, 2]-4) = 0, (m[2, 1]-1)*(m[2, 1]-2)*(m[2, 1]-3)*(m[2, 1]-4) = 0, (m[2, 2]-1)*(m[2, 2]-2)*(m[2, 2]-3)*(m[2, 2]-4) = 0, m[1, 1]*m[2, 2]-m[1, 2]*m[2, 1] = 0, (m[1, 1]-m[1, 2])*(m[1, 1]-m[2, 1])*(m[1, 1]-m[2, 2]) <> 0, (m[1, 2]-m[1, 1])*(m[1, 2]-m[2, 1])*(m[1, 2]-m[2, 2]) <> 0, (m[2, 1]-m[1, 1])*(m[2, 1]-m[1, 2])*(m[2, 1]-m[2, 2]) <> 0, (m[2, 2]-m[1, 1])*(m[2, 2]-m[1, 2])*(m[2, 2]-m[2, 1]) <> 0}

(3)

isolve(rels);  # no solution founds

 

A non systematic approach

Case of S(3)

 

 

M := Matrix(3, 3, symbol=m):
iM := {indices(M)}:

rels :=
  {
    det(M) = 0,

    # each term is equal to an integer between 1 and 9 included

    mul((M[1, 1]-k), k=1..9)=0,
    mul((M[1, 2]-k), k=1..9)=0,
    mul((M[1, 3]-k), k=1..9)=0,
    mul((M[2, 1]-k), k=1..9)=0,
    mul((M[2, 2]-k), k=1..9)=0,
    mul((M[2, 3]-k), k=1..9)=0,
    mul((M[3, 1]-k), k=1..9)=0,
    mul((M[3, 2]-k), k=1..9)=0,
    mul((M[3, 3]-k), k=1..9)=0,

    # the sum of all the terms is equal to 10*9/2 = 45

    M[1, 1]+M[1, 2]+M[1, 3]+M[2, 1]+M[2, 2]+M[2, 3]+M[3, 1]+M[3, 2]+M[3, 3]=45,

    # all the terms are different

    seq( mul(seq((M[op(im)]-M[op(ij)]), ij in iM minus {im})) <> 0, im in iM)
  }:

numelems(rels);

20

(4)

run_this := false:
if run_this then
  isolve(rels);  # requires a huge amount of memory and computational time
end if;

 


 

Download HowManyMatricesAreSingular.mw

Hi, 

I'm stucked in determining the intersection curve(s) of two (intersecting !) cylinders.
Plotting these curves can easily be done with plots:-intersectionplot, but I'm interested in finding the algebraic equations of this (these) curve(s).

I tried to do this while using either parametric or implicit representations of the two cylinders.

(For now on I'm using Maple 2015 and I wasn't capable to repoduce a few promising results I'd obtained at the office with Maple 2019 and parametric representations. So I mainly concentrated onimplicit representations).

If E(x,y,z) and E'(x,y,z) denote implicit representations of cylinders C and C', I had (naively) thought that simply solving 
E(x,y,z) = E'(x,y,z) with respect with x, y and z would have done the job.

Unfortunately, even for the simplest case of orthogonal circular cylinders of same radii, solve returns the couple of planes which contain the two intersection curves (ellipses) but not these ellipses themselves

Maybe there is a "simple" way to obtain the algebraic equation(s) of the intersection curve(s) but I wasn't capable to find it.
Instead of that I wrote a complicated stuff (please look to the attached file) which works well in some situations and not in others (see the last test case).

Could you please help me to answer this issue?

Thanks in advance

PS: no real need to consider tangent cylinders along a generatrix or one-point tangency. 

Intersecting_Cylinders.mw

Hi, 

While "numpoints=..." is said to be an admissible option of SurfaceOfRevolution (or VolumeOfRevolution, default value set to 50), the graphical results doesn't account for the value of numpoints.

Is this a known issue?

TIA

Surface0fRevolution.mw



 

Hi,

It seems that plottools:-extrude doesn't support the "style" option: no error returned, just the extusion being always of surface style.
Am I correct ?

1 2 3 4 5 6 7 Last Page 1 of 18