Advanced Search

vv

13992 Reputation

20 Badges

10 years, 38 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

M...

vv 13992
December 08 2017
1 0

To obtain directly the wanted result just use:

f:=10000/(1+30762*0.478^x)+5:
maximize(diff(f,x), location);

    1845.361367, {[{x = 14.00001597}, 1845.361367]}

Now let me explain why those complex numbers appeared.
The graph of the derivative is

Due to inherent round-offs (because you did not use rationals and solve instead of fsolve) your assigned M could be a bit bigger than the actual maximum of f'. So, the equation f'(x)=M will have only complex roots. To get rid of the imaginary part you may take Re(a), where a is the complex root.

Eigenvectors...

vv 13992
December 05 2017
1 2

Your matrix T is a standard test for eigenvector computation. It is symmetric but has “repeated” eigenvalues.

Maple can find the exact eigenvalues/eigenvectors in this case.

V,Q:=Eigenvectors(T);

But the test is designed for approximate computations (floating point).
So, use:
V1,Q1 := Eigenvectors(evalf(T));

(Then take Re to get rid of the small imaginary parts due to round-offs).

Note that in both cases, the matrix Q (or Q1)  containing the eigenvectors is not orthogonal (even if such a matrix exists for symmetric matrices). Eigenvectors() finds it only for float symmetric matrices with distinct eigenvalues.
If  orthogonality is needed then use GramSchmidt.
 

 

A simpler example to practice...

vv 13992
December 03 2017
1 3

I'd suggest to take a simpler example of a integro-diff system:

 

> 

e1:=diff(f(x),x) = f(x)*x - (x^2+1)*int(f(t),t=0..x)*g(x) + (1/2)*x^6+x^4-(1/2)*x^2+1;
e2:=diff(g(x),x) = f(x)*g(x)^2  -x^5-2*x^3+x;
ic:=f(1)=1, g(1)=2;

diff(f(x), x) = f(x)*x-(x^2+1)*(int(f(t), t = 0 .. x))*g(x)+(1/2)*x^6+x^4-(1/2)*x^2+1

  

diff(g(x), x) = f(x)*g(x)^2-x^5-2*x^3+x

  

f(1) = 1, g(1) = 2

 (1)
> 

dsolve({e1,e2,ic}, {f(x),g(x)}, numeric);

Error, (in dsolve/numeric/process_input) input system must be an ODE system, got independent variables {t, x}

  

Hence, not accepted. We transform it into an ODE system introducing F(x) = Int(f(t), t=0..x)  ==>

> 

E1:=diff(f(x),x) = f(x)*x-(x^2+1)*F(x)*g(x)+(1/2)*x^6+x^4-(1/2)*x^2+1;
E2:=diff(g(x),x) = f(x)*g(x)^2-x^5-2*x^3+x;
E3:=diff(F(x),x) = f(x);
IC:= f(1) = 1, g(1) = 2, F(0) = 0;

diff(f(x), x) = f(x)*x-(x^2+1)*F(x)*g(x)+(1/2)*x^6+x^4-(1/2)*x^2+1

  

diff(g(x), x) = f(x)*g(x)^2-x^5-2*x^3+x

  

diff(F(x), x) = f(x)

  

f(1) = 1, g(1) = 2, F(0) = 0

 (2)
> 

dsolve({E1,E2,E3,IC}, {f(x),g(x),F(x)}, numeric);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

  
> 

 

dsolve accepts it but is not able to solve using the implicit method.

 

You can play with it using different methods (already known by you from the other thread).

This system has the advantage that here we know an exact solution:  f(x) = x,  g(x) = x^2+1,  F(x)=x^2/2;

 

> 

simplify(eval([E1,E2,E3,IC], [f = (x->x),  g=(x -> x^2+1),  F=(x->x^2/2)]));

[1 = 1, 2*x = 2*x, x = x, 1 = 1, 2 = 2, 0 = 0]

 (3)
> 

 

Good look!


 

Download ode-ex.mw

rsolve...

vv 13992
December 03 2017
3 0

Denoting F(n) = f(0)+...+f(n) :

rsolve({f(n+1)=F(n), f(n)=F(n)-F(n-1), F(0)=1}, {f(n),F(n)});

(f(n) for n>0 of course)

First forty derivatives at 0...

vv 13992
December 01 2017
3 3 
p:=ln((1+x)*y) + exp((x^2)*(y^2)) - x - cos(x):
n:=40: Order:=n+1: interface(rtablesize=infinity):
Y:=convert(solve(series(p,x),y),polynom):
Vector[column](n, k ->  (D@@k)(y)(0)=eval(diff(Y,x$k),x=0));

(D(y))(0) = 0
((D@@2)(y))(0) = -2
((D@@3)(y))(0) = -2
((D@@4)(y))(0) = 55
((D@@5)(y))(0) = 96
((D@@6)(y))(0) = -3951
((D@@7)(y))(0) = -14986
((D@@8)(y))(0) = 599292
((D@@9)(y))(0) = 3756040
((D@@10)(y))(0) = -150583195
((D@@11)(y))(0) = -1451365992
((D@@12)(y))(0) = 56622924821
((D@@13)(y))(0) = 787829738580
((D@@14)(y))(0) = -29495018843850
((D@@15)(y))(0) = -570264628995262
((D@@16)(y))(0) = 20174284393035491
((D@@17)(y))(0) = 528406101043111104
((D@@18)(y))(0) = -17368560725263003371
((D@@19)(y))(0) = -607502621004862751222
((D@@20)(y))(0) = 18165803871525900350124
((D@@21)(y))(0) = 844845810439069562531252
((D@@22)(y))(0) = -22331822067408858287791607
((D@@23)(y))(0) = -1391102518495130466232499136
((D@@24)(y))(0) = 31137730503244680996639906697
((D@@25)(y))(0) = 2661297665528307975243093795336
((D@@26)(y))(0) = -46909537305923595629791853669346
((D@@27)(y))(0) = -5812207204576583472301320482485946
((D@@28)(y))(0) = 69478963661245363181701637264539327
((D@@29)(y))(0) = 14236970345056860768024681386350762560
((D@@30)(y))(0) = -71299313756801875315745344773429106215
((D@@31)(y))(0) = -38343398748183919497538363485168620149282
((D@@32)(y))(0) = -148767292026687176665730760219959487334644
((D@@33)(y))(0) = 110608691261039295372096825657213070145581600
((D@@34)(y))(0) = 1784186834273881762868359210688981525024869565
((D@@35)(y))(0) = -327307575749373732034049337639072936118525057752
((D@@36)(y))(0) = -11596412918978194063904912726445275723035187765539
((D@@37)(y))(0) = 899618765531255422817716130916175428672344561662044
((D@@38)(y))(0) = 67838929477109472349986391474357426166573238064107894
((D@@39)(y))(0) = -1466744386949476812393980478836845534302116051221857142
((D@@40)(y))(0) = -389148909149217184540519589123668475583410103807659331061

 

vertices...

vv 13992
December 01 2017
2 0
> 

planes:=2*sqrt(2)*x-2*z+sqrt(2), 2*sqrt(2)*x+2*z-sqrt(2), 2*sqrt(2)*y-2*z-sqrt(2), 2*sqrt(2)*y+2*z+sqrt(2);

2*2^(1/2)*x-2*z+2^(1/2), 2*2^(1/2)*x+2*z-2^(1/2), 2*2^(1/2)*y-2*z-2^(1/2), 2*2^(1/2)*y+2*z+2^(1/2)

 (1)
> 

vertices:=seq( eval([x,y,z], solve({planes[[({1,2,3,4}minus{i})[]]]}) ), i=1..4);

[1, 0, -(1/2)*2^(1/2)], [-1, 0, -(1/2)*2^(1/2)], [0, -1, (1/2)*2^(1/2)], [0, 1, (1/2)*2^(1/2)]

 (2)
> 

plots:-display(plottools:-tetrahedron([vertices]), transparency=0.5);

 
> 

 

 

 

parse...

vv 13992
November 30 2017
0 0

I prefer:

deq := FunctionAdvisor(DE, hypergeom([a, b], [c], z))[2, 1]:
Deq:=parse(convert(deq,string)):
dsolve(Deq, f(z));

 

select, diff...

vv 13992
November 29 2017
0 1 
T:=u-> degree( eval(u,diff=1), [f(x),g(x)]):
select( u -> T(u)<2, p);

 

Digits...

vv 13992
November 28 2017
0 2

In evalf[20], a = .1  is computed with Digits=10.
You should use subs(a=1/10) or Digits:=20

Explore...

vv 13992
November 26 2017
0 0 
# ex := 2.*ca*g+ka*ma+ka*meqc+keqc*ma-1.*ma*meqc*(ca*keqc+2.*g*ka)/(ca*ma+ca*meqc+g*ma): #TableauRouth[3,1]
ex := 
(2*ca^3*g*keqc*ma+2*ca^3*g*keqc*meqc+4*ca^2*g^2*ka*ma+4*ca^2*g^2*ka*meqc+2*ca^2*g^2*keqc*ma+ca^2*keqc^2*ma^2+4*ca*g^3*ka*ma+2*ca*g*ka^2*ma^2+4*ca*g*ka^2*ma*meqc+2*ca*g*ka^2*meqc^2+ca*g*ka*keqc*ma^2-3*ca*g*ka*keqc*ma*meqc+ca*g*keqc^2*ma^2+2*g^2*ka^2*ma^2-2*g^2*ka^2*ma*meqc+g^2*ka*keqc*ma^2)*(2*ca^2*g*ma+2*ca^2*g*meqc+2*ca*g^2*ma+ca*ka*ma^2+2*ca*ka*ma*meqc+ca*ka*meqc^2+ca*keqc*ma^2+g*ka*ma^2-g*ka*ma*meqc+g*keqc*ma^2): #[4,1]
Explore( g=solve(ex>0,g), parameters=[ca=0..5.0, ka=0..5.0, keqc=0..5.0, ma=0..5.0, meqc=0..5.0], echoexpression=false );

 

divergent...

vv 13992
November 26 2017
1 2

heunB is 1 for r=0 and continuous. The integral is divergent due to 1/r.

SemiAlgebraic...

vv 13992
November 26 2017
1 0

The inequalities are polynomial (after eliminating the denominators: a/b>0 <==> a*b>0).
So they can be solved with SemiAlgebraic or even solve, but a HUGE and useless result should be expected.

restart...

vv 13992
November 25 2017
0 2

Put restart in a separate execution group.

You may also switch to Assume using

with(MathematicalFunctions,Assume);
Assume(s>0);

 

Domains...

vv 13992
November 25 2017
0 2

For working with polynomials over GF, the simplest method is to use the Domains package .

This package does not seem to be very popular, and even some links in the documentation are broken., e.g. the link Domains/example should be
?examples,Domains.

 

If you are interested only in the field Fp i.e. GF(p,1) then Domains is not necessary because GF(p,1) is Zp and the operator mod is enough.

 

For example, to compute the gcd of two polynomials in Z7[x]:

 

> 

a := x^5+5*x^4+6*x+2;
b := x^5+3*x^3+5*x^2+2*x+5;

x^5+5*x^4+6*x+2

  

x^5+3*x^3+5*x^2+2*x+5

 (1)
> 

Gcd(a,b) mod 7;

x^2+1

 (2)

Note that the gcd of these polynomials in Z[x] is 1

> 

gcd(a,b);

1

 (3)

Now let's compute the same gcd using Domains:

> 

restart;
> 

with(Domains);

---------------------- Domains version 1.0 ---------------------
Initially defined domains are Z and Q the integers and rationals
Abbreviations, e.g. DUP for DenseUnivariatePolynomial, also made

  

[AbelianGroup, AbelianMonoid, AlgebraicExtension, Array, CommutativeRing, CramersRule, DenseExponentVector, DenseUnivariatePolynomial, DifferentialField, DistributedMultivariatePolynomial, EuclideanAlgorithm, EuclideanDomain, ExpandedNormalFormQuotientField, ExponentVector, FactoredNormalFormQuotientField, Field, Finite, FiniteField, FrobeniusNormalForm, GaloisField, Gaussian, GaussianInteger, GcdDomain, Generator, GrobnerBasis, Group, Integers, IntegralDomain, LazyUnivariatePowerSeries, List2Vector, MacaulayExponentVector, Maple, MapleDistributedMultivariatePolynomial, MapleExponentVector, MapleGaussianArithmetic, MapleMultivariatePolynomial, MapleUnivariatePolynomial, Matrix, MatrixDomain, Modp1Polynomial, Modp2Polynomial, Monoid, New, NormalForm, ODESolve, OrderedAbelianGroup, OrderedAbelianMonoid, OrderedDomain, OrderedField, OrderedSet, OrderedUnivariatePolynomial, PartiallyOrderedSet, PolynomialHomomorphism, Poset, PowerRemainder, PrimeExponentVector, PrincipalIdeal, Q, QuotientField, RandMatrix, RationalFunction, Rationals, RelativelyPrimeHeuristic, RepeatedSquaring, Ring, SemiGroup, Set, SingleExtensionGaloisField, Singular, Solve, SparseDistributedMultivariatePolynomial, SparseMultivariatePolynomial, SquareFreeNorm, SquareMatrix, TableDistributedMultivariatePolynomial, TranscendentalFunctions, UniqueFactorizationDomain, UnivariatePolynomial, UnivariatePowerSeries, Vector, Z, Zmod, addCategory, addProperties, addProperty, defOperation, defOperations, evaldomains, hasCategories, hasCategory, hasOperation, hasOperations, hasProperties, hasProperty, isDomain, isImplemented, newCategory, newDomain, notImplemented, show]

 (4)
> 

C := DenseUnivariatePolynomial(GF(7,1), x):
> 

a := x^5+5*x^4+6*x+2;

x^5+5*x^4+6*x+2

 (5)
> 

b := x^5+3*x^3+5*x^2+2*x+5;

x^5+3*x^3+5*x^2+2*x+5

 (6)
> 

a1:=C[Input](a);

`domains/DenseUnivariatePolynomial/badge0`(2, 6, 0, 0, 5, 1)

 (7)
> 

b1:=C[Input](b);

`domains/DenseUnivariatePolynomial/badge0`(5, 2, 5, 3, 0, 1)

 (8)
> 

C[Gcd](a1,b1);

`domains/DenseUnivariatePolynomial/badge0`(1, 0, 1)

 (9)

Note trhat the conversion to a1 and b1 is necessary:

> 

C[Gcd](a,b);

`domains/DenseUnivariatePolynomial/badge0`(1)

 (10)
> 

lprint(a);

x^5+5*x^4+6*x+2

  
> 

lprint(a1);

`domains/DenseUnivariatePolynomial/badge0`(2, 6, 0, 0, 5, 1)

  

 

A last remark. If you want to implement your own gcd, then note that the elements of GF(p,k) are themselves polynomials (actually in the quotient ring Zp[alpha]/(r), where r is an irreducible polynomial in the variable alpha, of degree k).  So, those  zppolys are elements of GF, not the polynomials in x].

 

By change of variables...

vv 13992
November 23 2017
0 7
> 

rel:=simplify(tan((x-t)/(alpha*(t^2+1)))-t);

-tan((-x+t)/(alpha*(t^2+1)))-t

 (1)
> 

alpha:=1/2;

1/2

 (2)
> 

X:=solve(rel,x);

(1/2)*t^2*arctan(t)+t+(1/2)*arctan(t)

 (3)
> 

t1:=fsolve(X=1);

.6177384029

 (4)
> 

evalf(Int(1/(1+t^2)*diff(X,t),t=0..t1));

.8903158470

 (5)
> 

 

Edit: without a change (i.e. directly):

restart;
rel:=simplify(tan((x-t)/(alpha*(t^2+1)))-t):
alpha:=1/2:
T:=solve(rel,t):
Int(1/(1+T^2),x=0..1):
evalf(%);

       .8903158471

First 80 81 82 83 84 85 86 Last Page 82 of 120