Experience 7 (Undetermined Coefficients with Sage).sagews

# 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)

$\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)
simplify(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 (note that since A0, B0 did not appear above they are zero)
solve ([4*B2 == -4, -4*A2 == -8, 2*A2 + 2*B1 == 8, -2*A1 + 2*B2 == -8], A1 ,B1, A2, B2)

[[A1 == 3, B1 == 2, 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)

$\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)

$\displaystyle x \ {\mapsto}\ -x$

$\displaystyle x \ {\mapsto}\ -{\left(x + 1\right)} e^{\left(-x\right)}$

$\displaystyle x \ {\mapsto}\ -x^{2} - x - 1$