SharedExperience 7 (Undetermined Coefficients with Sage).sagewsOpen in CoCalc
# Method of Undetermined Coefficients
# Example y ’’-2y ’+y=-3-x+x^2
x=var("x")
y= function ("y")(x)
# For RHS of ODE = -3-x+x^2 we try quadratic yp
A=var("A")
B=var("B")
C=var("C")
yp=A+B*x+C*x^2
expand(diff(yp,x,2) - 2*diff(yp,x) + yp) # LHS of ODE y’’-2y’+y  - Equate the results of this with the RHS
C*x^2 + B*x - 4*C*x + A - 2*B + 2*C
# Solves system to find coefficients
solve ([C==1 , B -4*C==-1, A -2*B +2*C== -3] ,A,B,C)
[[A == 1, B == 3, C == 1]]

# restates particular solution and checks it
yp =1+3* x+1*x^2
show (yp)
expand(diff(yp,x,2) - 2*diff(yp,x) + yp)
x2+3x+1\displaystyle x^{2} + 3 \, x + 1
x^2 - x - 3

# Example y ’’+y=(8 -4*x)*cos (x) -(8 -8*x)*sin(x)
x=var("x")
y= function ("y")(x)
# For RHS of ODE = (8 -4*x)*cos(x) -(8 -8*x)*sin(x) 
# Normally you would just use trig functions times quadratics 
# but since the general solution is already of this form we need to use a higher order polynomial
# to get a linearly independent particular solution
A0= var("A0")
B0= var("B0")
A1= var("A1")
B1= var("B1")
A2= var("A2")
B2= var("B2")
yp =(A0+A1*x+A2*x^2)*cos(x)+(B0+B1*x+B2*x^2)*sin(x)
expand(diff(yp,x,2) + yp) #LHS of ODE y ’’ + y Equate the results of this with RHS
4*B2*x*cos(x) - 4*A2*x*sin(x) + 2*A2*cos(x) + 2*B1*cos(x) - 2*A1*sin(x) + 2*B2*sin(x)
# Solves system to find coefficients
solve ([4*B2 == -4, -4*A2 == 8, 2*A2 + 2*B1 == 8, -2*A1 + 2*B2 == -8] ,A1 ,B1, A2, B2)
[[A1 == 3, B1 == 6, A2 == -2, B2 == -1]]
# restates particular solution and checks it
yp = (3*x - 2*x^2)*cos(x) + (6*x - x^2)*sin(x)
show (yp)
expand(diff(yp,x,2) + yp)
(2x23x)cos(x)(x26x)sin(x)\displaystyle -{\left(2 \, x^{2} - 3 \, x\right)} \cos\left(x\right) - {\left(x^{2} - 6 \, x\right)} \sin\left(x\right)
-4*x*cos(x) + 8*x*sin(x) + 8*cos(x) - 8*sin(x)
# Variation of Parameters
# Example (x -1)y ’’-xy ’+y=(x -1) ^2 where y1=x, y2=e^x solve homogeneous ODE
# we try for a solution of the form yp= u1y1 + u2y2
u1p =var("u1p")
u2p =var("u2p")
# We require that u1p *y1+u2p*y2 =0 and u1p *y1 ’+ u2p*y2 ’=F/Po
solve ([u1p*x + u2p*e^x == 0 , u1p + u2p*e^x == x-1], u1p ,u2p )
[[u1p == -1, u2p == x*e^(-x)]]
u1p (x)=-1 # Define u1 ’ as a function
u1(x)=u1p. integrate (x) # Integrate to get u1
show (u1)
u2p (x)=x*e^(-x) # define u2 ’ as a function
u2(x)=u2p. integrate (x) # Integrate to get u2
show (u2)
# Show particular solution
yp= simplify ( expand (u1*x+u2*e^x))
show (yp)
x  x\displaystyle x \ {\mapsto}\ -x
x  (x+1)e(x)\displaystyle x \ {\mapsto}\ -{\left(x + 1\right)} e^{\left(-x\right)}
x  x2x1\displaystyle x \ {\mapsto}\ -x^{2} - x - 1