Mariusz Iwaniuk

1571 Reputation

14 Badges

9 years, 289 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by Mariusz Iwaniuk

You made a  simple mistake in the equation (18).You  wrote xi(1),but  it should be: 1/xi(1).

See attached file:

feijisuan.mw

 

You forgot to deliver a inital conditions to DEplot function.

with(DEtools):

DE3 := {diff(y(x), x) = y(x)-z(x), diff(z(x), x) = z(x)-2*y(x)};

DEplot(DE3, [y(x), z(x)], x = 0 .. 3, y = 0 .. 2, z = -4 .. 4, [[y(0) = 1, z(0) = 1]], arrows = medium, numpoints = 200);

 

Download DEplot.mw

For more information execute: ?DEplot in  Maple command line.

On Maple 2018.1  definition doublefactorial is the same as 2016 nothing has changed.

From Maple help: ?factorial;

"Note: In Maple, !! is used for repeated factorials and so it does not indicate the double factorial."

then use: doublefactorial(n);

Maybe you can use alias function:

alias(`bi-factorial` = doublefactorial);

`bi-factorial`(15);

#2027025

doublefactorial(15);

#2027025

 

What version of Maple do you use?.

If you are using an older version, the latest version of Maple int correctly solves.


 

NULL

restart

infolevel[IntegrationTools] := 3

lambda[1] := 3/10; lambda[2] := 2*(1/10); alpha := 4

3/10

 

1/5

 

4

(1)

int(2*alpha^2*Z*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity)

Definite Integration:   Integrating expression on Z=0..infinity
Definite Integration:   Using the integrators [distribution, piecewise, series, o, polynomial, ln, lookup, cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper, asymptotic, ftoc, ftocms, meijerg, contour]

IntegralTransform LookUp Integrator:   Integral might be a laplace transform with expr=32/(exp(1/5*Z)+3)*Z/(exp(3/10*Z)+3)^2 and s=-3/10.

Definite Integration:   Integrating expression on T=0..infinity
Definite Integration:   Using the integrators cook
Cook LookUp Integrator:   but does not fit into the first class
Cook LookUp Integrator:   returning answer from cook pattern 1a
Definite Integration:   Returning integral unevaluated.
Definite Integration:   Integrating expression on T=0..infinity
Definite Integration:   Using the integrators cook

Cook LookUp Integrator:   but does not fit into the first class
Cook LookUp Integrator:   returning answer from cook pattern 1a
Definite Integration:   Returning integral unevaluated.
Definite Integration:   Integrating expression on T=0..infinity
Definite Integration:   Using the integrators cook
Definite Integration:   Returning integral unevaluated.

IntegralTransform LookUp Integrator:   Integral transform returned unevaluated.
LookUp Integrator:   unable to find the specified integral in the table

Definite Integration:   Method ftoc succeeded.
Definite Integration:   Finished sucessfully.

 

(1600/27)*ln(2)-(400/9)*dilog(1+(1/3)*3^(2/3))-(1000/81)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(1/3)-(1000/81)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(1/3)+(800/81)*3^(5/6)*arctan((2/3)*3^(1/6)-(1/3)*3^(1/2))-(400/81)*3^(5/6)*Pi-(400/27)*3^(1/6)*Pi+(250/81)*Pi^2+(250/27)*ln(3)^2-(400/81)*3^(1/3)*ln(3^(2/3)-3^(1/3)+1)+(800/27)*3^(1/6)*arctan((2/3)*3^(1/6)-(1/3)*3^(1/2))-((200/27)*I)*3^(1/6)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))+((200/27)*I)*3^(1/6)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-((1000/81)*I)*3^(5/6)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))+((1000/81)*I)*3^(5/6)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-((200/9)*I)*3^(1/2)*dilog(1-((1/3)*I)*3^(1/2))+((200/9)*I)*3^(1/2)*dilog(1+((1/3)*I)*3^(1/2))+(800/81)*3^(1/3)*ln(1+3^(1/3))+(200/81)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(2/3)+(2000/81)*dilog(1+(1/3)*3^(2/3))*3^(1/3)-(400/81)*dilog(1+(1/3)*3^(2/3))*3^(2/3)+(400/81)*ln(3^(2/3)-3^(1/3)+1)*3^(2/3)-(800/81)*ln(1+3^(1/3))*3^(2/3)+(200/81)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(2/3)+(200/3)*dilog(1-((1/3)*I)*3^(1/2))+(200/3)*dilog(1+((1/3)*I)*3^(1/2))-(400/9)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-(400/9)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-(4000/729)*3^(5/6)*Pi*ln(3)-(800/243)*3^(1/6)*Pi*ln(3)-(800/729)*3^(2/3)*Pi^2+(4000/729)*3^(1/3)*Pi^2+(100/9)*3^(1/2)*ln(3)*Pi

(2)

evalf(((200/27)*I)*3^(1/6)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))+((1000/81)*I)*3^(5/6)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))+((200/9)*I)*3^(1/2)*dilog(1+((1/3)*I)*3^(1/2))-(4000/729)*3^(5/6)*Pi*ln(3)-(800/243)*3^(1/6)*Pi*ln(3)+(100/9)*3^(1/2)*ln(3)*Pi-((200/27)*I)*3^(1/6)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-((1000/81)*I)*3^(5/6)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-((200/9)*I)*3^(1/2)*dilog(1-((1/3)*I)*3^(1/2))+(1600/27)*ln(2)-(400/9)*dilog(1+(1/3)*3^(2/3))-(400/9)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))-(400/9)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))+(250/81)*Pi^2+(250/27)*ln(3)^2+(200/3)*dilog(1-((1/3)*I)*3^(1/2))+(200/3)*dilog(1+((1/3)*I)*3^(1/2))-(400/81)*3^(1/3)*ln(3^(2/3)-3^(1/3)+1)+(800/27)*3^(1/6)*arctan((2/3)*3^(1/6)-(1/3)*3^(1/2))+(800/81)*3^(1/3)*ln(1+3^(1/3))+(200/81)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(2/3)+(2000/81)*dilog(1+(1/3)*3^(2/3))*3^(1/3)-(400/81)*dilog(1+(1/3)*3^(2/3))*3^(2/3)+(400/81)*ln(3^(2/3)-3^(1/3)+1)*3^(2/3)-(800/81)*ln(1+3^(1/3))*3^(2/3)+(200/81)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(2/3)-(800/729)*3^(2/3)*Pi^2+(4000/729)*3^(1/3)*Pi^2-(1000/81)*dilog(1-((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(1/3)-(1000/81)*dilog(1+((1/2)*I)*3^(1/6)-(1/6)*3^(2/3))*3^(1/3)+(800/81)*3^(5/6)*arctan((2/3)*3^(1/6)-(1/3)*3^(1/2))-(400/81)*3^(5/6)*Pi-(400/27)*3^(1/6)*Pi)``

19.31389182+0.*I

(3)

``


Maybe you can try and see what hapens ?

infolevel[IntegrationTools] := 3;
lambda[1] := 3/10; lambda[2] := 2*(1/10); alpha := 4:
int(2*alpha^2*Z*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity, method = ftocms);

Download aquestion_ver_3.mw

 

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

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