11838 Reputation

7 years, 186 days

not constant...

It is not constant as

plots[odeplot](solA0,[theta,r(theta)],theta=varepsilon..5/1000);

shows. It is only a precision problem, abserr and relerr are used probably with some heuristic tests which are not enough when the solution is singular.

@tomleslie  Is seems to be hard to ...

Is seems to be hard to obtain o prescribed precision.

In your worksheet one obtains the value of A:

A0 := 199.961318896644828500644760320

and the solution

Now, if we increase a bit the precision in dsolve without parameters for this A0:

solA0 := dsolve( {Ode1, r(varepsilon) = A0, (D(r))(varepsilon) = 0},
numeric,
range = varepsilon .. 5*(1/1000),
abserr=1.0e-14,
relerr=1.0e-13
);

one obtains for solA0(5/1000);

i.e.  r(theta) has only 8 correct digits!

How many digits can we guarantee for A0?

@tomleslie I do not think that the optim...

@tomleslie

I do not think that the optimization is essential here; fsolve will do it after adjusting Digits.

The real problem is the accuracy of dsolve with parameters.

precision?...

I think we should have serious concerns about the precision because of the singularity.

For example fsolve (instead of NLPSolve) refuses to solve the equation even in the interval  199.96..199.97.

If we set Digits:=27 before fsolve ==>

199.961306129991433460330154

but the problem is how good is the approximation given by dsolve with parameters.

@acer  I did not do it on purpose. ...

I did not do it on purpose.

showstat(plottools:-`transform/object`);

simply does not include it.

Anyway I don't think that maplesoft will have something to lose, on the contrary!

@taizoon  So all you know about the...

So all you know about theta(z,p)  is contained in    pde=0.
lambda(p) seems to be a given function.
The substitution p=0 is of no use.

You may consider theta = theta(z) as depending of the parameter p.

So,

ode := diff(theta(z),z,z) + lambda(p)^2*theta(z)+p*lambda(p)^2*(sin(theta(z))-theta(z));
s:=dsolve(ode=0,theta(z));

# it contains 2 solutions (given implicitely). Take e.g. the first one.

s[1];

keep in mind that the constants _C1 and _C2 depend on p.

@taizoon  You did not post the whol...

You did not post the whole problem and I cannot guess what you are trying to do.

@Markiyan Hirnyk  The simple ideea ...

The simple ideea is to use:

min{ |f(x)| + |g(x)| : x in X } =
min{ a + b : -a <= f(x) <= a,  -b <= g(x) <= b,  x in X, a in R, b in R}.

So, the standard simplex can be used if f,g are affine
(and provides an  _exact_ solution).

@Carl Love MultiInt:=proc(A,B)local...

MultiInt:=proc(A,B)
local r:=NULL,x,Ax,Bx;
for x in {op(A)} intersect {op(B)} do
Ax:=select(`in`,A,[x]);   # or     Ax:=select(u->evalb(u=x), A)
Bx:=select(`in`,B,[x]);
if nops(Ax)<nops(Bx) then r:=r,op(Ax) else r:=r,op(Bx) fi
od;
[r]
end;

Probably the efficiency is worse (it would be interesting to test) but you must admit that it is much easier to read/understand.

common...

It does not work correctly
L([1,2,3,1,1,1,2],[1,1,2,2,2,2,2,2]) ;
should return [1,1,2,2]

Probably a "one line definition" for L is more subtle. Of course, a longger one is easy.

lists...

If multipicity is considered in the lists, one should use e.g.

L:=(A,B) -> select(u->member(u,B), A);

(actually this should be modified for the correct number of common multiplicity)

Thank you for the info. I know that Geog...

Thank you for the info. I know that Geogebra is also popular. They include CASes, but it would be nice take advantage of a full featured system such as Maple.

It is not dificult to program it.Here is...

It is not dificult to program it.

Here is a not optimal one where the circles are represented as lists [x,y,r].

myapo:=proc(c1,c2,c3)
local i1,i2,i3,x,y,r,s,u,c,apo:=NULL;
for i1 in [-1,1] do for i2 in [-1,1] do for i3 in [-1,1] do
s:=solve({ (x-c1[1])^2+(y-c1[2])^2 = (r*i1 + c1[3])^2,
(x-c2[1])^2+(y-c2[2])^2 = (r*i2 + c2[3])^2,
(x-c3[1])^2+(y-c3[2])^2 = (r*i3 + c3[3])^2 },
{x,y,r});
for u in [s] do
c:=evalf([eval(x,u),eval(y,u),eval(r,u)]);
if not(type(c,list(realcons))) or (c[3]<0) then next fi;
apo:=apo,c od;
od od od;
print(cat("number of circles:",nops([apo])));
apo
end:

To test it with your data:

c1,c2,c3:=[0,0,5],[5,4,2],[13,0,3]:

c8:=myapo(c1,c2,c3);
ccc:=[seq(circle([u[1],u[2]],u[3],color=red),u=[c1,c2,c3]), seq(circle([u[1],u[2]],u[3],color=black),u=[c8])]:
with(plottools):with(plots):
display(ccc);

try-catch...

You may use

try
A := Apollonius( c1, c2, c3 )
catch:   A:= MyApolonius(c1,c2,c3)
end try;

which catches any exception.

 First 158 159 160 161 162 Page 160 of 162
﻿