using SymPy x = Sym("x") @syms a b c a,b,c = Sym("a,b,c") u = symbols("u") x = symbols("x", real=true) y1, y2 = symbols("y1, y2", positive=true) alpha = symbols("alpha", integer=true, positive=true) solve(x^2 + 1) # ±i are not real solve(y1 + 1) # -1 is not positive @syms u1 positive=true u2 positive=true solve(u1 + u2) # empty, though solving u1 - u2 is not. PI, E, oo [asin(1), asin(Sym(1))] @syms x y ex = x^2 + 2x + 1 subs(ex, x, y) subs(ex, x, 0) x,y,z = symbols("x,y,z") ex = x + y + z subs(ex, (x,1), (y,pi)) subs(ex, x=>1, y=>pi) ex(x=>1, y=>pi) ex(1, pi) ex |> subs(x, 1) ex |> replace(y, pi) N(PI) # floating-point value evalf(PI) N(PI, 30) evalf(PI, 30) @syms x y z p = x^2 + 3x + 2 factor(p) expand(prod([(x-i) for i in 1:5])) factor(x^2 - 2) q = x*y + x*y^2 + x^2*y + x collect(q, x) collect(q, y) simplify(q) expand((x + 1)*(x - 2) - (x - 1)*x) factor(exp(2x) + 3exp(x) + 2) r = 1/x + 1/x^2 together(r) apart( (4x^3 + 21x^2 + 10x + 12) / (x^4 + 5x^3 + 5x^2 + 4x)) top = (x-1)*(x-2)*(x-3) bottom = (x-1)*(x-4) top/bottom r = expand(top) / expand(bottom) cancel(r) @syms x y nonnegative=true a real=true simplify(x^a * y^a - (x*y)^a) x,y,a = symbols("x,y,a") simplify(x^a * y^a - (x*y)^a) powsimp(x^a * y^a - (x*y)^a, force=true) theta = symbols("theta", real=true) p = cos(theta)^2 + sin(theta)^2 simplify(p) simplify(sin(2theta) - 2sin(theta)*cos(theta)) expand_trig(sin(2theta)) a,b,c,x = symbols("a, b, c, x") p = a*x^2 + b*x + c coeff(p, x^2) # a coeff(p, x) # b p(x=>0) Sym[[coeff(p, x^i) for i in N(degree(p)):-1:1]; p(x=>0)] coeffs(p) # fails q = Poly(p, x) coeffs(q) real_roots(x^2 - 2) p = (x-3)^2*(x-2)*(x-1)*x*(x+1)*(x^2 + x + 1) real_roots(p) solve(p) polyroots(p) p = x^5 - x + 1 polyroots(p) rts = solve(p) [N(r) for r in rts] # or map(N, rts) nroots(p) solve(cos(x) - sin(x)) u = solveset(cos(x) - sin(x)) v = solveset(x^2 - 4) elements(v) solve(cos(x) - x) nsolve(cos(x) - x, 1) a,b,c = symbols("a,b,c", real=true) p = a*x^2 + b*x + c solve(p, x) solveset(p, x) solve(p) x, y = symbols("x,y", real=true) exs = [2x+3y-6, 3x-4y-12] d = solve(exs) map(ex -> subs(ex, d), exs) a,b,c,h = symbols("a,b,c,h", real=true) p = a*x^2 + b*x + c fn = cos exs = [fn(0*h)-p(x=>0), fn(h)-p(x => h), fn(2h)-p(x => 2h)] d = solve(exs, [a,b,c]) quad_approx = subs(p, d) n = 3 x, c = symbols("x,c") as = Sym["a$i" for i in 0:(n-1)] bs = Sym["b$i" for i in 0:(n-1)] p = sum([as[i+1]*x^i for i in 0:(n-1)]) q = sum([bs[i+1]*(x-c)^i for i in 0:(n-1)]) solve(p-q, bs) solve(Eq(x, 1)) solve(x ⩵ 1) x, y = symbols("x,y", real=true) exs = [2x+3y ⩵ 6, 3x-4y ⩵ 12] ## Using \Equal[tab] d = solve(exs) x = symbols("x") using Plots # plot(x^2 - 2, -2,2) plot(sin(2x), cos(3x), 0, 4pi) limit(sin(x)/x, x, 0) limit(sin(x)/x, x=>0) limit((1+1/x)^x, x => oo) a = symbols("a", positive=true) ex = (sqrt(2a^3*x-x^4) - a*(a^2*x)^(1//3)) / (a - (a*x^3)^(1//4)) ex(x=>a) # or subs(ex, x, a) subs(denom(ex), x, a), subs(numer(ex), x, a) limit(ex, x => a) quad_approx limit(quad_approx, h => 0) limit(sign(x), x => 0) limit(sign(x), x => 0, dir="-"), limit(sign(x), x => 0, dir="+") f(x) = sin(5x)/x limit(f, 0) f(x) = 1/x^(log(log(log(log(1/x)))) - 1) hs = [10.0^(-i) for i in 6:16] ys = [f(h) for h in hs] [hs ys] limit(f(x), x, 0, dir="+") x, h = symbols("x,h") f(x) = exp(x)*sin(x) limit((f(x+h) - f(x)) / h, h, 0) diff(f(x), x) diff(x^x, x) diff(exp(-x^2), x, 2) diff(exp(-x^2), x, x, x) # same as diff(..., x, 3) f(x) = (12x^2 - 1) / (x^3) diff(f(x), x) |> solve f(x) = exp(x)*cos(x) diff(f, 2) x,y = symbols("x,y") ex = x^2*cos(y) Sym[diff(ex,v1, v2) for v1 in [x,y], v2 in [x,y]] ex = Derivative(exp(x*y), x, y, 2) doit(ex) F, G = SymFunction("F"), SymFunction("G") F(x) diff(F(x)) x,y = symbols("x, y") ex = y^4 - x^4 - y^2 + 2x^2 ex1 = ex(y=>F(x)) ex2 = diff(ex1, x) ex3 = solve(ex2, F'(x))[1] ex4 = ex3(F(x) => y) w, h, P = symbols("w, h, P", nonnegative=true) r = w/2 A = w*h + 1//2 * (pi * r^2) p = w + 2h + pi*r h0 = solve(P-p, h)[1] A1 = A(h => h0) coeffs(Poly(A1, w)) coeff(collect(expand(A1), w), w^2) A1(w => 0) b = solve(subs(P-p, h, 0), w)[1] c = solve(diff(A1, w), w)[1] atb = A1(w => b) atc = A1(w => c) atc - atb simplify(atc - atb) integrate(x^3, x) n = symbols("n", real=true) ex = integrate(x^n, x) ex(n => 3) integrate(x^2, (x, 0, 1)) integrate(x^5*sin(x), x) f(x) = (x^4 + x^2 * exp(x) - x^2 - 2x*exp(x) - 2x - exp(x)) * exp(x) / ( (x-1)^2 * (x+1)^2 * (exp(x) + 1) ) f(x) integrate(f(x), x) x, y = symbols("x,y") integrate(x*y, (y, 0, 1), (x, 0, 1)) integrate(x^2*y, (y, 0, sqrt(1 - x^2)), (x, -1, 1)) integ = Integral(sin(x^2), x) doit(integ) f(x) = exp(x) * cos(x) integrate(f) integrate(sin, 0, pi) s1 = series(exp(sin(x)), x, 0, 4) s2 = series(cos(exp(x)), x, 0, 6) simplify(s1 * s2) removeO(s1) i, n = symbols("i, n") summation(i^2, (i, 1, n)) sn = Sum(1/i^2, (i, 1, n)) doit(sn) limit(doit(sn), n, oo) x,y = symbols("x,y") v = [1,2,x] w = [1,y,3] dot(v,w) cross(v,w) ex = x^2*y - x*y^2 Sym[diff(ex,var) for var in [x,y]] Sym[diff(ex, v1, v2) for v1 in [x,y], v2 in [x,y]] hessian(ex) x,y = symbols("x,y") M = [1 x; x 1] diagm(ones(Sym, 5)) M^2 det(M) eigvecs(M) M[:eigenvects]() M = Sym[1 2 3 0 0; 4 10 0 0 1] vs = nullspace(M) [M*vs[i] for i in 1:3] M = [1 x; x 1] P, D = diagonalize(M) # M = PDP^-1 D, M - P*D*inv(P) F = SymFunction("F") diffeq = Eq(diff(F(x), x, 2) - 2*diff(F(x)) + F(x), sin(x)) ex = dsolve(diffeq) diffeq = F''(x) - 2F'(x) + F(x) - sin(x) dsolve(diffeq) ex1 = rhs(ex) solve(ex1(x => 0), Sym("C1")) ex2 = ex1(Sym("C1") => -1//2) solve( subs(diff(ex2, x), x, 0) - 1, Sym("C2") ) ex3 = ex2(Sym("C2") => 3//2) t, m,k,alpha = symbols("t,m,k,alpha") v = SymFunction("v") ex = Eq( (m/k)*v'(t), alpha^2 - v(t)^2 ) ex[:classify_ode]() dsolve(ex) y = SymFunction("y") a, x = symbols("a,x") eqn = y'(x) - 3*x*y(x) - 1 x0, y0 = 0, 4 out = dsolve(eqn, x, (y, x0, y0)) u = rhs(out) diff(u, x) - 3*x*u - 1 x0, y0 = 0, a out = dsolve(eqn, x, (y, x0, y0)) as = -2:0.6:2 ex = rhs(out) p = plot(ex(a=>as[1]), -1.8, 1.8, ylims=(-4, 4)) for i in as[2:end] plot!(p, ex(a=>i), -1.8, 1.8, ylims=(-4, 4)) end p y = SymFunction("y") x = symbols("x") eqn = y''(x) + 5y'(x) + 6y(x) out = dsolve(eqn, x, (y, 0, 1), (y', 0, 1)) plot(rhs(out), -1/3, 2) eqn = y''(x) + y(x) - exp(x) dsolve(eqn, x, (y, 0, 1), (y, 1, 1//2))