Mariusz Iwaniuk

1526 Reputation

14 Badges

9 years, 123 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by Mariusz Iwaniuk

 

Substitution x = t^2:

IntegrationTools:-Change(Int(1/(1+sqrt(x)), x = 0 .. 1), x = t^2);

value(%);

#2-2*ln(2)

Or using convert(expr,elementary):

Int(1/(1+sqrt(x)), x = 0 .. 1) = int(1/(1+sqrt(x)), x = 0 .. 1);

#Int(1/(1+sqrt(x)), x = 0 .. 1) = MeijerG([[0, 0, 1/2], []], [[1/2, 0], [-1]], 1)/Pi

Substitution: 1=x

convert(MeijerG([[0, 0, 1/2], []], [[1/2, 0], [-1]], x)/Pi, elementary);

#-ln(1-x)/x-2*arctanh(sqrt(x))/x+2/sqrt(x)

limit(%, x = 1);

#2-2*ln(2)

It's a bug in int !


 

restart; f := ((1/2-I*t)^(-s)-(1/2+I*t)^(-s))/(2*I); fc := evalc(f); f1 := `assuming`([simplify(int(f, t = 0 .. infinity))], [s > 1]); f2 := `assuming`([simplify(int(fc, t = 0 .. infinity))], [s > 1])

-((1/2)*I)*((1/2-I*t)^(-s)-(1/2+I*t)^(-s))

 

exp(-(1/2)*s*ln(1/4+t^2))*sin(s*arctan(2*t))

 

2^(s-1)/(s-1)

 

4^s/(2*s-2)

(1)

eval(f1, s = 3)

2

(2)

int(eval(fc, s = 3), t = 0 .. infinity, numeric)

2.

(3)

eval(f2, s = 2)

8

(4)


 

Download function_evaluation_goes_wrong_v2.mw

 

Only for: z>0,j>0,Re(z),Re(j)


 

restart

func := expand((-exp(j*z)*Ei(1, j*z)+I*Pi*exp(-j*z)+exp(-j*z)*Ei(1, -j*z))/(2*j))

-(1/2)*exp(j*z)*Ei(1, j*z)/j+((1/2)*I)*Pi/(j*exp(j*z))+(1/2)*Ei(1, -j*z)/(j*exp(j*z))

(1)

func2 := eval(func, [Ei(1, j*z) = -Ei(-j*z), Ei(1, -j*z) = -Ei(j*z)])

(1/2)*exp(j*z)*Ei(-j*z)/j+((1/2)*I)*Pi/(j*exp(j*z))-(1/2)*Ei(j*z)/(j*exp(j*z))

(2)

plot([Re(eval(func, [j = 1/2])), Re(eval(func2, [j = 1/2]))], z = 0 .. 5, color = ["red", "black"])

 

eval(Ei(1, -z) = -Ei(z)+(1/2)*ln(z)-(1/2)*ln(1/z)-ln(-z), z = j*z)

Ei(1, -j*z) = -Ei(j*z)+(1/2)*ln(j*z)-(1/2)*ln(1/(j*z))-ln(-j*z)

(3)

Ei(z) = -Ei(1, -z)+(1/2)*ln(z)-(1/2)*ln(1/z)-ln(-z)

func3 := simplify(eval(func, [Ei(1, j*z) = -Ei(-j*z)+(1/2)*ln(-j*z)-(1/2)*ln(-1/(j*z))-ln(j*z), Ei(1, -j*z) = -Ei(j*z)+(1/2)*ln(j*z)-(1/2)*ln(1/(j*z))-ln(-j*z)]))

(1/4)*(exp(j*z)*ln(-1/(j*z))-ln(1/(j*z))*exp(-j*z)+((2*I)*Pi-2*ln(-j*z)-2*Ei(j*z)+ln(j*z))*exp(-j*z)-exp(j*z)*(ln(-j*z)-2*ln(j*z)-2*Ei(-j*z)))/j

(4)

plot([Re(eval(func, [j = 1/2])), Re(eval(func3, [j = 1/2]))], z = 0 .. 5, color = ["red", "black"])

 

help("Ei # you can find formula here")

``


 

Download Ei.mw

interface(typesetting = extended);
F := n-> ifactors(ithprime(n)-2)[2, 1, 1] ;
F(11);
print(F);

#                               29
#n -> ifactors(ithprime(n) - 2)[2, 1, 1]

 

On Maple 2018.1

I a little speed up computation.The plot is now gererated by Maple about 1 minute.

If N=0 then integral is probably singular ,I changed q = 0.1e-2 to q = 0.1e-1.


 

restart; Omega := 2*Pi*N; R0 := a*tanh((a^2-mu)/(2*T_c))*ln((2*a^2+2*a*q+q^2-2*mu-I*Omega)/(2*a^2-2*a*q+q^2-2*mu-I*Omega))/q-2; T_c := 0.169064e-1; mu := .869262; N := 10; R1 := unapply(Int(R0, a = 0.1e-3 .. 100, method = _Gquad), q)

2*Pi*N

 

a*tanh((1/2)*(a^2-mu)/T_c)*ln((2*a^2+2*a*q+q^2-2*mu-(2*I)*Pi*N)/(2*a^2-2*a*q+q^2-2*mu-(2*I)*Pi*N))/q-2

 

0.169064e-1

 

.869262

 

10

 

proc (q) options operator, arrow; Int(a*tanh(29.57459897*a^2-25.70807505)*ln((2*a^2+2*a*q+q^2-1.738524-(20*I)*Pi)/(2*a^2-2*a*q+q^2-1.738524-(20*I)*Pi))/q-2, a = 0.1e-3 .. 100, method = _Gquad) end proc

(1)

Digits := 5

n := 5; plots:-pointplot({seq([q, evalf(abs(R1(q)))], q = 0.1e-1 .. 10, 1/n)}, connect = true)

 

``


 

Download PLOT_fun.mw

A simple example than yours and  we convert to piecewise to understand more:

convert(signum(0, x), piecewise);

#piecewise(x < 0, -1, x = 0, 0, 0 < x, 1)

 

 

convert(signum(0, x,1), piecewise);

#piecewise(x < 0, -1, 0 <= x, 1)

 

Using derivatives:

diff(signum(0, x),x);

convert(%, piecewise);

#signum(1, x)

#piecewise(x = 0, undefined, 0)

 

 

diff(signum(0, x, 1),x);

convert(%, piecewise);

#signum(1, x,1)

#piecewise(x = 0, undefined, 0)

 

_Envsignum0 := 0; convert(signum(0, x), piecewise);

#piecewise(x < 0, -1, x = 0, 0, 0 < x, 1)

 

_Envsignum0 := 1; convert(signum(0, x), piecewise);

#piecewise(x < 0, -1, 0 <= x, 1)

 

_Envsignum0 := 2; convert(signum(0, x), piecewise);

#piecewise(x < 0, -1, x = 0, 2, 0 < x, 1)

 

_Envsignum0 := -2; convert(signum(0, x), piecewise);

#piecewise(x < 0, -1, x = 0, -2, 0 < x, 1)

 

1.

convert(Ei(1, -ln(a-2)), Sum, dummy = n);

#-gamma-ln(-ln(a-2))+Sum((-1)^n*(-ln(a-2))^(n+1)/(factorial(n+1)*(n+1)), n = 0 .. infinity)

2.

convert(Ei(a, b), Sum, dummy = n);

eval(%, b = 1);

value(%);

limit(%, a = 1);

evalf(%);

#0.2193839344

 

Probably your equation have no analytical solution.


 

restart

-3*Pi^2*3^(2/3)*4^(1/3)*S^3*(epsilon/Pi)^(2/3)+16*a*(S/Pi)^(3/2)*Pi^5+24*S^4*Pi = 0

-3*Pi^2*3^(2/3)*4^(1/3)*S^3*(epsilon/Pi)^(2/3)+16*a*(S/Pi)^(3/2)*Pi^5+24*S^4*Pi = 0

(1)

simplify(-3*Pi^2*3^(2/3)*4^(1/3)*S^3*(epsilon/Pi)^(2/3)+16*a*(S/Pi)^(3/2)*Pi^5+24*S^4*Pi = 0)

-3*3^(2/3)*2^(2/3)*S^3*Pi^(4/3)*epsilon^(2/3)+16*a*Pi^(7/2)*S^(3/2)+24*S^4*Pi = 0

(2)

solve((2),[S]);

[[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5)^2]]

(3)

allvalues([[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5)^2]])

[[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5, index = 1)^2]], [[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5, index = 2)^2]], [[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5, index = 3)^2]], [[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5, index = 4)^2]], [[S = 0], [S = RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*epsilon^(2/3)*_Z^3+16*Pi^(7/2)*a+24*Pi*_Z^5, index = 5)^2]]

(4)

epsilon := 1; a := 1

1

 

1

(5)

rhs((3)[2][1])

RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*_Z^3+16*Pi^(7/2)+24*Pi*_Z^5)^2

(6)

[evalf(allvalues(RootOf(-3*3^(2/3)*2^(2/3)*Pi^(4/3)*_Z^3+16*Pi^(7/2)+24*Pi*_Z^5)^2))]

[HFloat(1.071956098421639)+HFloat(2.5245757717629815)*I, -HFloat(1.932416441934298)-HFloat(1.5609688301792792)*I, HFloat(2.9299146528586517), -HFloat(1.932416441934298)+HFloat(1.5609688301792792)*I, HFloat(1.071956098421639)-HFloat(2.5245757717629815)*I]

(7)

``


 

Download rootof_ver2.mw

 


 

restart

eq:=I*mu*A(t[2])*omega[0]*(1/2)-A(t[2])*omega[0]*(diff(B(t[2]), t[2]))+I*(diff(A(t[2]), t[2]))*omega[0]-(1/4)*A(t[2])^5*beta[2]*omega[0]^2-(1/4)*A(t[2])^3*beta[1]*omega[0]^2-(1/2)*F[0]*exp(I*sigma*t[2]-I*B(t[2]))+5*alpha[2]*A(t[2])^5*(1/16)+3*alpha[1]*A(t[2])^3*(1/8);

((1/2)*I)*mu*A(t[2])*omega[0]-A(t[2])*omega[0]*(diff(B(t[2]), t[2]))+I*(diff(A(t[2]), t[2]))*omega[0]-(1/4)*A(t[2])^5*beta[2]*omega[0]^2-(1/4)*A(t[2])^3*beta[1]*omega[0]^2-(1/2)*F[0]*exp(I*sigma*t[2]-I*B(t[2]))+(5/16)*alpha[2]*A(t[2])^5+(3/8)*alpha[1]*A(t[2])^3

(1)

eq2:=simplify(eq);

-(1/2)*F[0]*exp(-I*(-sigma*t[2]+B(t[2])))+I*(diff(A(t[2]), t[2]))*omega[0]+(1/2)*(-2*(diff(B(t[2]), t[2]))*omega[0]+(-(1/2)*beta[2]*omega[0]^2+(5/8)*alpha[2])*A(t[2])^4+(-(1/2)*beta[1]*omega[0]^2+(3/4)*alpha[1])*A(t[2])^2+I*mu*omega[0])*A(t[2])

(2)

subs(-sigma*t[2]+B(t[2]) = C(t2), eq2)

-(1/2)*F[0]*exp(-I*C(t2))+I*(diff(A(t[2]), t[2]))*omega[0]+(1/2)*(-2*(diff(B(t[2]), t[2]))*omega[0]+(-(1/2)*beta[2]*omega[0]^2+(5/8)*alpha[2])*A(t[2])^4+(-(1/2)*beta[1]*omega[0]^2+(3/4)*alpha[1])*A(t[2])^2+I*mu*omega[0])*A(t[2])

(3)

eval(eq2,-sigma*t[2]+B(t[2])=C(t2));

-(1/2)*F[0]*exp(-I*C(t2))+I*(diff(A(t[2]), t[2]))*omega[0]+(1/2)*(-2*(diff(B(t[2]), t[2]))*omega[0]+(-(1/2)*beta[2]*omega[0]^2+(5/8)*alpha[2])*A(t[2])^4+(-(1/2)*beta[1]*omega[0]^2+(3/4)*alpha[1])*A(t[2])^2+I*mu*omega[0])*A(t[2])

(4)

 


 

Download no_change_ver2.mw

The standard way to solve systems of equations numerically is using fsolve.  For example:

fsolve({x^2-y-3, x*y+2*y^2-4 = 0}, {x = 0.5 .. 2.5, y = 0 .. 3});

#{x = 2.000000000, y = 1.000000000}

 

That's what the fsolve command will do with it.

infolevel[fsolve]:=6:
fsolve({x^2-y-3, x*y+2*y^2-4 = 0}, {x = 0.5 .. 2.5, y = 0 .. 3});


#fsolve: trying multivariate Newton iteration
#fsolve:
#guess vector Vector(2, [2.1564724756843477,.13446305023563622])
#fsolve: norm of errors: 5.1897839975614897
#fsolve: new norm: 3.6816464474611895
#fsolve: iter = 1 |incr| = 1.4002 new values x = 2.1216 y = 1.4998
#fsolve: new norm: .43227630491510698
#fsolve: iter = 2 |incr| = .53690 new values x = 2.0189 y = 1.0655
#fsolve: new norm: 0.9718475514730430e-2
#fsolve: iter = 3 |incr| = 0.82477e-1 new values x = 2.0005 y = 1.0015
#fsolve: new norm: 0.5297608520287e-5
#fsolve: iter = 4 |incr| = 0.19416e-2 new values x = 2.0000 y = 1.0000
#fsolve: new norm: 0.1568000e-11
#fsolve: iter = 5 |incr| = 0.10595e-5 new values x = 2.0000 y = 1.0000
#fsolve: new norm: 0.
#fsolve: iter = 6 |incr| = 0.31360e-12 new values x = 2.0000 y = 1.0000
#               {x = 2.000000000, y = 1.000000000}

Or using code from:https://www.mapleprimes.com/questions/200377-Nonlinear-Equations-With-Newoton-Newton#answer201481 thanks for Carl Love.


 

restart:

F:= unapply(<x1^2-y1-3,x1*y1+2*y1^2-4>, x1, y1):#Equations

NewtonsMethod:= proc(F, x0, {maxiters::posint:= 99}, {epsilon::positive:= 10^(3-Digits)})
local
     x, y,
     J:= unapply(VectorCalculus:-Jacobian(F(x,y), [x,y]), x, y),
     X:= x0,
     newX:= <1,1>,
     err:= <1+epsilon, epsilon>,
     k
;
     for k to maxiters while LinearAlgebra:-Norm(err)/LinearAlgebra:-Norm(newX) > epsilon do
          newX:= X - LinearAlgebra:-LinearSolve(J(X[1],X[2]), F(X[1],X[2]));
          err:= newX - X;
          X:= newX
     end do;
     if k > maxiters then  WARNING("Did not converge.")  end if;
     X
end proc:

NewtonsMethod(F, <0.5, 3.0>);

Vector(2, {(1) = 2.000000000000003, (2) = .9999999999999993})

(1)


 

Download MultiNewton.mw

 

 

 


 

restart

K := 1; k[e] := 1; d := 1; h := 1

F := proc (s, x) options operator, arrow; -2*K*(s+K)^2*k[e]*sinh((s+K)*x/d)/(s*d^2*(4*(s+K)^3*cosh((s+K)*h/d)*h/d^2-4*K*(s+K)*cosh((s+K)*h/d)*h/d+4*sinh((s+K)*h/d)*K)) end proc

F(s, x)

-2*(s+1)^2*sinh((s+1)*x)/(s*(4*(s+1)^3*cosh(s+1)-4*(s+1)*cosh(s+1)+4*sinh(s+1)))

(1)

NULL

Digits := 14

NLaplace := proc (F, t, x, n) local h, theta, z, dz; h := 2*Pi/n; theta := -Pi+(k+1/2)*h; z := n*(.5017*theta*cot(.6407*theta)-.6122+.2645*I*theta)/t; dz := n*((-1)*.5017*.6407*theta*csc(.6407*theta)^2+.5017*cot(.6407*theta)+.2645*I)/t; Re(evalf(-((1/2)*I)*h*(sum(exp(z*t)*F(z, x)*dz, k = 0 .. n))/Pi)) end proc

invLAP := proc (t, x) options operator, arrow; NLaplace(F, t, x, 32) end proc

plot(invLAP(1/2, x), x = -1 .. 1)

 

plot3d(invLAP(t, x), x = -1 .. 1, t = 0 .. 2)

 

NULL


 

Download ILT_(1).mw

This package is build-in Maple,you don't need to download anything:

Execute:

with(Student[Basics]);

LinearSolveSteps((x+1)/(2*y*z) = 4*y^2/z+3*x/y, x);

LinearSolveSteps function have some limitation only  for linear equations give steps.

LinearSolveSteps("x^2 -2*x+1 = 0", x); #Nonlinear equation. It doesn't work.

 

 

Eliminating variable m1 and m2 from equations, then we have only eq. with 2 variables k1 and k2.

Ploting we see they do not intersect anywhere, even at point {0,0}.


 

restart

sys := {(k1*m1+k1*m2+k2*m1+sqrt(k1^2*m1^2+2*k1^2*m1*m2+k1^2*m2^2+2*k1*k2*m1^2-2*k1*k2*m1*m2+k2^2*m1^2))/(2*m1*m2) = (2*Pi*8.78)^2, .8347192842*k1/m1 = (2*Pi*4.8515)^2, -(-30.85287127*k1-k2)/m2 = (2*Pi*8.78)^2, -(-k1*m1-k1*m2-k2*m1+sqrt(k1^2*m1^2+2*k1^2*m1*m2+k1^2*m2^2+2*k1*k2*m1^2-2*k1*k2*m1*m2+k2^2*m1^2))/(2*m1*m2) = (2*Pi*4.8515)^2}

{.8347192842*k1/m1 = 929.2055780, -(-30.85287127*k1-k2)/m2 = 3043.328048, -(1/2)*(-k1*m1-k1*m2-k2*m1+(k1^2*m1^2+2*k1^2*m1*m2+k1^2*m2^2+2*k1*k2*m1^2-2*k1*k2*m1*m2+k2^2*m1^2)^(1/2))/(m1*m2) = 929.2055780, (1/2)*(k1*m1+k1*m2+k2*m1+(k1^2*m1^2+2*k1^2*m1*m2+k1^2*m2^2+2*k1*k2*m1^2-2*k1*k2*m1*m2+k2^2*m1^2)^(1/2))/(m1*m2) = 3043.328048}

(1)

eq2 := eliminate(sys, {m1, m2})

[{m1 = 0.8983149735e-3*k1, m2 = 0.1013787235e-1*k1+0.3285876462e-3*k2}, {-0.2219750257e-1*k1^2-0.2848636636e-3*k1*k2+(1/2)*(0.1217974307e-3*k1^4-0.9347355846e-5*k1^3*k2+0.3245892274e-6*k1^2*k2^2)^(1/2), 0.2944183889e-2*k1^2-0.3391728651e-3*k1*k2+(1/2)*(0.1217974307e-3*k1^4-0.9347355846e-5*k1^3*k2+0.3245892274e-6*k1^2*k2^2)^(1/2)}]

(2)

plots:-implicitplot([eq2[2, 1] = 0, eq2[2, 2] = 0], k1 = -1 .. 1, k2 = -10 .. 10, view = [-1 .. 1, -10 .. 10], gridrefine = 5, color = ["Red", "Blue"], rational, thickness = 3)

 

plots:-implicitplot([eq2[2, 1] = 0, eq2[2, 2] = 0], k1 = -0.1e-3 .. 0.1e-3, k2 = -0.1e-1 .. 0.1e-1, view = [-0.1e-3 .. 0.1e-3, -0.1e-1 .. 0.1e-1], gridrefine = 5, color = ["Red", "Blue"], rational, thickness = 3)

 

``


 

Download No_solutions.mw

You must add 4 more (intial,boundary) condition if you what to solve by numeric method.
e.g. :{ U(0, t) = 1, U(1, t) = 0, V(0, t) = 0, V(1, t) = 0}


 

PDESYS := [diff(U(x, t), t)-(diff(U(x, t), x, x))-2*U(x, t)*(diff(U(x, t), x))+diff(U(x, t)*V(x, t), x), diff(V(x, t), t)-(diff(V(x, t), x, x))-2*V(x, t)*(diff(V(x, t), x))+diff(U(x, t)*V(x, t), x)]

[diff(U(x, t), t)-(diff(diff(U(x, t), x), x))-2*U(x, t)*(diff(U(x, t), x))+(diff(U(x, t), x))*V(x, t)+U(x, t)*(diff(V(x, t), x)), diff(V(x, t), t)-(diff(diff(V(x, t), x), x))-2*V(x, t)*(diff(V(x, t), x))+(diff(U(x, t), x))*V(x, t)+U(x, t)*(diff(V(x, t), x))]

(1)

NULL

IBC := [U(x, 0) = sin(x), V(x, 0) = sin(x), U(0, t) = 1, U(1, t) = 0, V(0, t) = 0, V(1, t) = 0]

[U(x, 0) = sin(x), V(x, 0) = sin(x), U(0, t) = 1, U(1, t) = 0, V(0, t) = 0, V(1, t) = 0]

(2)

sol := pdsolve(PDESYS, IBC, numeric)

_m489011961280

(3)

p1 := sol:-plot(V(x, t), t = 1/8, numpoints = 100, x = 0 .. 1, color = ["Blue"], legend = ["V(x,t)"])

p2 := sol:-plot(U(x, t), t = 1/20, numpoints = 100, x = 0 .. 1, color = ["Green"], legend = ["U(x,t)"])

plots:-display({p1, p2})

 

sol:-plot3d(U(x, t), t = 0 .. 1, x = 0 .. 1, grid = [100, 100])

 

sol:-plot3d(V(x, t), t = 0 .. 1, x = 0 .. 1, grid = [100, 100])

 

``


Executed in Maple 2018 !!!

Download pdesys.mw

First 11 12 13 14 15 16 17 Page 13 of 20