A power series can be written as \( y = \displaystyle \sum_{n=0}^{\infty} c_n x^n \).
The first derivative would be \( y' = \displaystyle \sum_{n=1}^{\infty} n c_n x^{n-1} \). Note that this starts at \( n = 1 \) since the \( n = 0 \) term was a constant and disappeared when you took the derivative.
The second derivative would be \( y'' = \displaystyle \sum_{n=2}^{\infty} n (n-1) c_n x^{n-2} \)
Maxima uses brackets [] for subscripts, so \( c_n \) would be written as c[n]. The maxima code for \( \infty \) is inf. That allows us to write the three values above as:
Y:sum(c[n]*x^n,n,0,inf)$
Y1:sum(n*c[n]*x^(n-1),n,1,inf)$
Y2:sum(n*(n-1)*c[n]*x^(n-2),n,2,inf)$Let's say that we wanted to solve the differential equation \( (x+4) y'' -x y' + 3y = \sin(2x) \)
After entering the definitions above, here are the next set of steps for working out the problem in Maxima.
(x+4)*Y2-x*Y1+3*Y$
expand(%)$
intosum(%)$
factorsum(%);(x+4)*Y2-x*Y1+3*Y enters the left hand side of the differential equation. We will hold off on the right hand side for the time being.expand(%) distributes the (x+4) and creates two separate summations.intosum(%) moves the x, 4, -x, and 3 into the summations.factorsum(%) simplifies the insides of the summations to combine any powers of x that might exist. You could also use expand() or factor(), but in this case, factor() factored a negative sign out of the entire expression while factorsum() didn't.This is the output from Maxima.
\[-\left( \sum_{n=1}^{\infty }n\,{c}_{n}\,{x}^{n}\right) +3\,\left( \sum_{n=0}^{\infty }{c}_{n}\,{x}^{n}\right) +\left( \sum_{n=2}^{\infty }\left( n-1\right) \,n\,{c}_{n}\,{x}^{n-1}\right) +4\,\sum_{n=2}^{\infty }\left( n-1\right) \,n\,{c}_{n}\,{x}^{n-2}\]
The first step in combining summations is to get everything so that it has a common power on the x. We write this as \( x^k \).
Let's say that we wanted to do this by hand for the term \( 4\,\sum_{n=2}^{\infty }\left( n-1\right) \,n\,{c}_{n}\,{x}^{n-2}\).
Right now, the power on x is \( n-2 \) and we would like it to be \( k \). So we use the substitution \( k = n-2\) or the equivalent \( n = k + 2\)
Make those substitutions into the sum:
\[4\,\sum_{\color{blue}{k+2}=2}^{\infty }\left( \color{blue} {k+2}-1\right) \,(\color{blue} {k+2})\,{c}_{\color{blue} {k+2}}\,{x}^{\color{blue} {k}}\]
and that simplifies to be
\[ 4\, \sum_{k=0}^{\infty} (k+1)(k+2)c_{k+2}x^k \]
There is a Maxima command to do all that for us. It's changevar() and takes four arguments.
k-n+2.So, if you entered this:
changevar(4*sum((n-1)*n*c[n]*x^(n-2),n,2,inf),k-n+2,k,n);Then you would get this:
\[4\,\sum_{k=0}^{\infty }\left( {k}^{2}+3\,k+2\right) \,{c}_{k+2}\,{x}^{k}\]
And that is the same thing we got, just with the \((k+1)(k+2)\) multiplied out to get \( k^2+3k+2\)..
The part() command will give you, among other things, the individual terms of an expression. We can use it to grab the fourth part of the last answer by using this expression: part(%,4).
That means that we could have written the changevar() command as:
changevar(part(%,4),k-n+2,k,n);Note that this should be part of the same block as the four steps above, right after the factorsum(), or % will refer to the wrong thing.
It is absolutely crucial that you look at the terms before you start using the part() command as Maxima will reorder them as it does calculations.
Here is the expression we had before, but now shown with the substitutions that we need to make:
\[\underbrace{-\left( \sum_{n=1}^{\infty }n\,{c}_{n}\,{x}^{n}\right) +3\,\left( \sum_{n=0}^{\infty }{c}_{n}\,{x}^{n}\right)}_{\color{blue}{k=n}} +\underbrace{\left( \sum_{n=2}^{\infty }\left( n-1\right) \,n\,{c}_{n}\,{x}^{n-1}\right)}_{\color{blue}{k=n-1}} +\underbrace{4\,\sum_{n=2}^{\infty }\left( n-1\right) \,n\,{c}_{n}\,{x}^{n-2}}_{\color{blue}{k=n-2}}\]
Here is the command needed to do this. Again, make sure this is immediately after the factorsum() command from before.
changevar(part(%,1)+part(%,2),k-n,k,n)+
changevar(part(%,3),k-n+1,k,n)+
changevar(part(%,4),k-n+2,k,n);The resulting expression is
\[4\,\left( \sum_{k=0}^{\infty }\left( {k}^{2}+3\,k+2\right) \,{c}_{k+2}\,{x}^{k}\right) +\left( \sum_{k=1}^{\infty }\left( {k}^{2}+k\right) \,{c}_{k+1}\,{x}^{k}\right) -\left( \sum_{k=1}^{\infty }k\,{c}_{k}\,{x}^{k}\right) +3\,\sum_{k=0}^{\infty }{c}_{k}\,{x}^{k}\]
Notice how all of the summations have the same \( x^k \) in them.
The problem now is that some summations begin with \( k = 0 \) and some begin with \( k = 1 \). We need them to start at the same place.
The trick to doing this is to pull the \( k = 0 \) terms out of those summations so they start with \( k = 1 \) and then leave it as a summation.
This is using the identity that \( \displaystyle \sum_{k=0}^{\infty} a_k = a_0 + \sum_{k=1}^{\infty} a_k \).
That identity could be extended if needed. For example, if you needed to start with \( k = 3 \), then you could pull out three terms to get \( \displaystyle \sum_{k=0}^{\infty} a_k = a_0 + a_1 + a_2 + \sum_{k=3}^{\infty} a_k \).
In the first summation, \( 4\,\left( \sum_{k=0}^{\infty }\left( {k}^{2}+3\,k+2\right) \,{c}_{k+2}\,{x}^{k}\right) \), when you find the \( k = 0 \) term, you get \( 4 (0^2+3(0)+2))c_{0+2} x^0 = 8 c_2 \)
In the fourth summation, \( 3\,\sum_{k=0}^{\infty }{c}_{k}\,{x}^{k}\), when you find the \( k = 0 \) term, you get \( 3 c_0 x^0 = 3 c_0 \).
\[8c_2 + 4\,\left( \sum_{k=1}^{\infty }\left( {k}^{2}+3\,k+2\right) \,{c}_{k+2}\,{x}^{k}\right) +\left( \sum_{k=1}^{\infty }\left( {k}^{2}+k\right) \,{c}_{k+1}\,{x}^{k}\right) -\left( \sum_{k=1}^{\infty }k\,{c}_{k}\,{x}^{k}\right) +3c_0 + 3\,\sum_{k=1}^{\infty }{c}_{k}\,{x}^{k}\]
Now all the summations start with \( k = 1 \) and the all have \( x^k \) in them, so we are ready to combine them.
List the terms you had to pull out separately and then combine the summations into a single summation.
\[ \underbrace{8c_2 + 3c_0}_{\color{blue}{k=0}} + \sum_{k=1}^{\infty} \left [ 4(k^2+3k+2)c_{k+2} + (k^2+k)c_{k+1}-kc_k+3c_k \right ] x^k \]
Simplifying (the last two terms inside the summation are both \( c_k \):
\[ 8c_2 + 3c_0 + \sum_{k=1}^{\infty} \left [ 4(k^2+3k+2)c_{k+2} + (k^2+k)c_{k+1}-(k-3)c_k \right ] x^k \]
Maxima will do these last two steps for you using the sumcontract() command. However, you must first bring all the coefficeints back inside the summations and you will probably want to factor it when you're done.
intosum(%)$
sumcontract(%);
LHS:factorsum(%);The result is the same thing we got by hand:
\[\left( \sum_{k=1}^{\infty }\left( \left( 4\,{k}^{2}+12\,k+8\right) \,{c}_{k+2}+\left( {k}^{2}+k\right) \,{c}_{k+1}+\left( 3-k\right) \,{c}_{k}\right) \,{x}^{k}\right) +8\,{c}_{2}+3\,{c}_{0}\]
The LHS label is so we can use that later.
The right hand side was \( \sin (2x) \), which can be written as a power series in the form of a Maclaurin series.
The Maclaurin series is \( \sin x = x - \frac {x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \dots \), so we make a subsitution to get \( \sin 2x = 2x - \frac {(2x)^3}{3!} + \frac{(2x)^5}{5!} - \frac{(2x)^7}{7!} + \dots \) or \( \sin 2x = 2x - \frac {4}{3}x^3 + \frac{4}{15}x^5 - \frac{8}{315}x^7 + \dots \).
Maxima will do this quickly for us:
RHS:taylor(sin(2*x),x,0,7);Undetermined Coefficients is a method we used with partial fractions. It involves setting similar terms from the left side equal to the rights side.
\( k = 0 \) corresponds to the constant terms and gives \( 8c_2 + 3c _0 = 0 \). Solving for the higher indexed term gives \( c_2 = -\frac{3}{8} c_0 \).
From here on, we're using the summation to get the terms. The coefficient on \( x^k \) is \( 4(k^2+3k+2)c_{k+2} + (k^2+k)c_{k+1}-(k-3)c_k \).
Maxima can do this, too! (whew!), but once again, you need to be careful with the part() command.
Make sure this works before you put it into block, it may change from problem to problem.
The summation is the first term, but we don't want the entire summation, we want the first part inside the summation, so we use:
REC:coeff(part(LHS,1,1),x,k);This assigns the coefficient on \( x^k \), which is \(\left( 4\,{k}^{2}+12\,k+8\right) \,{c}_{k+2}+\left( {k}^{2}+k\right) \,{c}_{k+1}+\left( 3-k\right) \,{c}_{k}\) to the label REC.
Unfortunately, that won't return the constant by asking for the coefficient on \(x^0\) because the summation doesn't have an \(x\) on the outside and gets included. Instead, we'll manually collect the constants with the part() command. Again, this is highly dependent on your problem!
CON:part(LHS,2)+part(LHS,3);Now, let's put it all together at the end of the factorsum() from the first part.
RHS:taylor(sin(2*x),x,0,7)$
REC:coeff(part(LHS,1,1),x,k);
CON:part(LHS,2)+part(LHS,3);
makelist(subst(k=n,REC)=coeff(RHS,x,n),n,1,6);
append(%,[CON=coeff(RHS,x,0)]);
SOL:linsolve(%,makelist(c[k],k,2,9));Here are what the lines accomplish:
Note, that if we had to pull out more than one term in the original summation, we could have appended more than just the constant.
We have two unknowns, \( c_0 \) and \( c_1 \). This means that we will be able to find two solutions, one in terms of \( c_0 \) and the second in terms of \( c_1 \). If you only have one unknown, then you only need to make one subsitution.
The first solution will involve \(c_0\) only. To do this, we substitute \(c_0 = 1\) and \(c_1 = 0\) into the solution.
subst([c[0]=1,c[1]=0],SOL);This gives us a list
\[{c}_{2}=-\frac{3}{8},{c}_{3}=\frac{11}{96},{c}_{4}=-\frac{5}{768},{c}_{5}=-\frac{241}{15360},{c}_{6}=\frac{59}{23040},{c}_{7}=\frac{1217}{1290240},{c}_{8}=-\frac{327}{2293760},{c}_{9}=-\frac{3911}{82575360}\]
These are the coefficients and the powers on the x's are the same as the subscript on the term. Don't forget that we had \(c_0 = 1\) and \( c_1 = 0\) as well, which means that there is a \( 1 + 0x = 1\) to add in.
Putting those back into a summation, we get the first solution
\[ y_1 = [1-\frac{3\,{x}^{2}}{8}+\frac{11\,{x}^{3}}{96}-\frac{5\,{x}^{4}}{768}-\frac{241\,{x}^{5}}{15360}+\frac{59\,{x}^{6}}{23040}+\frac{1217\,{x}^{7}}{1290240}-\frac{327\,{x}^{8}}{2293760}-\frac{3911\,{x}^{9}}{82575360}+\dots\]
The second solution will involve \(c_10\) only. To do this, we substitute \(c_0 = 0\) and \(c_1 = 1\) into the solution.
subst([c[0]=0,c[1]=1],SOL);This gives us a list
\[{c}_{2}=0,{c}_{3}=-\frac{1}{12},{c}_{4}=\frac{1}{96},{c}_{5}=-\frac{1}{640},{c}_{6}=\frac{1}{2880},{c}_{7}=-\frac{13}{161280},{c}_{8}=\frac{17}{860160},{c}_{9}=-\frac{461}{92897280}\]
These are the coefficients and the powers on the x's are the same as the subscript on the term. Don't forget that we had \(c_0 = 0\) and \( c_1 = 1\) as well, which means that there is a \( 0 + 1x = x\) to add in.
Putting those back into a summation, we get the second solution
\[y_2 = x-\frac{{x}^{3}}{12}+\frac{{x}^{4}}{96}-\frac{{x}^{5}}{640}+\frac{{x}^{6}}{2880}-\frac{13\,{x}^{7}}{161280}+\frac{17\,{x}^{8}}{860160}-\frac{461\,{x}^{9}}{92897280} + \dots\]
The general solution to the differential equation is the linear combination \( y = c_1 y_1 + c_2 y_2 \), where \(c_1\) and \( c_2 \) turn are arbitrary constants that we can't find unless we have a specific point that hte curve passes through.
If you only had one solution, then the general solution would be \( y = c_1 y_1 \).
Okay, here is the complete Maxima code. The powerdisp:false$ tells Maxima to list the summations in order of decreasing powers and powerdisp:true$ tells it to list it in ascending power. We normally write polynomials in descending power, but power series in ascending power. This definitely affects the part() command. The default is false and that's how I worked the problem, but I changed it when displaying the solutions. Then, when I went back, it wouldn't work anymore because the parts weren't the same. I reset it to the default at the top and the bottom just so it didn't throw anything else off and then everything worked, even if you have to run it again.
You can copy and paste this to duplicate this example problem, but do NOT just change the beginning and expect it to work. The part() commands require knowing what position expressions are in. There are also a lot of $ at the end of lines that could be ; to display. Again, this is the end result, not the first time through.
powerdisp:false$
Y:sum(c[n]*x^n,n,0,inf)$
Y1:sum(n*c[n]*x^(n-1),n,1,inf)$
Y2:sum(n*(n-1)*c[n]*x^(n-2),n,2,inf)$
powerdisp:false$
(x+4)*Y2-x*Y1+3*Y$
expand(%)$
intosum(%)$
factorsum(%)$
changevar(part(%,1)+part(%,2),k-n,k,n)+
changevar(part(%,3),k-n+1,k,n)+
changevar(part(%,4),k-n+2,k,n)$
intosum(%)$
sumcontract(%)$
LHS:factorsum(%)$
RHS:taylor(sin(2*x),x,0,7)$
REC:coeff(part(LHS,1,1),x,k);
CON:part(LHS,2)+part(LHS,3)$
makelist(subst(k=n,REC)=coeff(RHS,x,n),n,1,7);
append(%,[CON=0])$
SOL:linsolve(%,makelist(c[k],k,2,9));
powerdisp:true$
subst([c[0]=1,c[1]=0],SOL)$
1+subst(%,sum(c[k]*x^k,k,2,9));
subst([c[0]=0,c[1]=1],SOL)$
x+subst(%,sum(c[k]*x^k,k,2,9));
powerdisp:false$