<h1>
    M 391 Differential Equations Final Project
</h1>
<h2>
    Due 11/18/20, 5 pm
</h2>
<h3>
Problem 1   
</h3>
By first setting $a_0=1$ and $a_1=0$ and then setting $a_0=0$ and $a_1=1$, find a pair of independent solutions to equations (1), (2), and (3) for p=0,1,2,3. $\newline$
Hermite's Equation (1) $$y''-2xy'+2py=0$$

Suppose $y = \sum_{n  = 0}^\infty a_n x^n$ solves (1) Then we have $y' = \sum_{n  = 1}^\infty na_n x^{n-1}$ and $y'' = \sum_{n  = 2}^\infty (n-1)na_n x^{n-2}$.  Plugging these into (1), we get:
\begin{align}
\sum_{n  = 2}^\infty (n-1)na_n x^{n-2} - 2x\sum_{n  = 1}^\infty na_n x^{n-1} + 2p \sum_{n  = 0}^\infty a_n x^n = 0  \\
\sum_{n  = 0}^\infty (n+2)(n+1)a_{n+2} x^{n} - \sum_{n  = 1}^\infty 2na_n x^{n} + \sum_{n  = 0}^\infty 2pa_n x^n = 0
\end{align}
Since the $n=0$ of the $y'$ term is $2*0*a_0*x^0 = 0$, we can write:
\begin{align}
\sum_{n  = 0}^\infty (n+2)(n+1)a_{n+2} x^{n} - \sum_{n  = 0}^\infty 2na_n x^{n} + \sum_{n  = 0}^\infty 2pa_n x^n = 0 \\
\sum_{n  = 0}^\infty ((n+2)(n+1)a_{n+2} - 2na_n + 2pa_n) x^n = 0 \\
(n+2)(n+1)a_{n+2} - 2na_n + 2pa_n = 0, \forall n \in \mathbb{N}
\end{align}
This gives us the recurrence relation:
\begin{align}
a_{n+2} = \frac{2(n-p)}{(n+2)(n+1)}a_n
\end{align}
If we set $a_0=1$ and $a_1=0$, we get $a_1 = a_3 = 0$. Thus we get the following coefficients:
\begin{align}
a_2 = \frac{2(0-p)}{2} \\
a_4 = \frac{2^2(2-p)(0-p)}{4!}\\
\end{align}
These initial conditions give us the independent solution $$g_p = a_0(1 + \frac{2(0-p)}{2}x^2 + \frac{2^2(2-p)(0-p)}{4!}x^4 + ...)$$
For $g_p$ where $p \in \{0,1,2,3\}$, we get the following solutions (see Problem 4 for more details):
\begin{align}
g_0 = 1\\
g_1 = 1 - x^2\\
g_2 = -(1/9)x^10 - (1/7)x^8 - (1/5)x^6 - (1/3)x^4 - x^2 + 1\\
g_3 = (2/7)x^10 + (3/7)x^8 + (4/5)x^6 + 3x^4 - 6x^2 + 1
\end{align}
If we set $a_0=0$ and $a_1=1$, we get $a_0 = a_2 = 0$. Thus, for $p \in \{0,1,2,3\}$ we get the following coefficients:
\begin{align}
a_3 = \frac{2(1-p)}{3!}\\
a_5 = \frac{2^2(3-p)(1-p)}{5!}
\end{align}
These initial conditions give us the independent solution $$h_p = a_1(x + \frac{2(1-p)}{3!}x^3 + \frac{2^2(3-p)(1-p)}{5!}x^5 + ...)$$
For $h_p$ where $p \in \{0,1,2,3\}$, we get the following solutions (see Problem 4 for more details):
\begin{align}
h_0 = (1/11)x^11 + (1/9)x^9 + (1/7)x^7 + (1/5)x^5 + (1/3)x^3 + x\\
h_1 = x\\
h_2 = -(2/33)x^11 - (5/63)x^9 - (4/35)x^7 - (1/5)x^5 - (2/3)x^3 + x\\
h_3 = -(2/3)x^3 + x
\end{align}

Legendre's equation (2) $$(1−x^2)y′′−2xy′+p(p+1)y=0$$
Suppose $y = \sum_{n  = 0}^\infty a_n x^n$ solves (2) Then we have $y' = \sum_{n  = 1}^\infty na_n x^{n-1}$ and $y'' = \sum_{n  = 2}^\infty (n-1)na_n x^{n-2}$.  Plugging these into (1), we get:
\begin{align}
\sum_{n  = 2}^\infty (n-1)na_n x^{n-2} - x^2\sum_{n  = 2}^\infty (n-1)na_n x^{n-2} - 2x\sum_{n  = 1}^\infty na_n x^{n-1} + p(p+1)\sum_{n=0}^\infty a_nx^n = 0 \\
\sum_{n = 0}^\infty (n+1)(n+2)a_{n+2}x^n - \sum_{n=0}^\infty (n-1)na_nx^n - \sum_{n  = 1}^\infty 2na_n x^n + \sum_{n=0}^\infty p(p+1)a_nx^n = 0 \\
\sum_{n=0}^\infty ((n+1)(n+2)a_{n+2} - (n-1)na_n - 2na_n + p(p+1)a_n)x^n = 0 \\
(n+1)(n+2)a_{n+2} - (n-1)na_n - 2na_n + p(p+1)a_n = 0
\end{align}
This gives the recurrence relation: $$a_{n+2} = \frac{n(n+1)-p(p+1)}{(n+1)(n+2)}$$

If we set $a_0=1$ and $a_1=0$, we get $a_1 = a_3 = 0$. Thus we get the following coefficients:
\begin{align}
a_2 = -\frac{p(p+1)}{2} \\
a_4 = \frac{p(p-2)(p+3)(p+1)}{4!}\\
\end{align}
These initial conditions give us the independent solution $$g_p = a_0(1 - \frac{p(p+1)}{2}x^2 + \frac{p(p-2)(p+3)(p+1)}{4!}x^4 - ...) = a_0 \sum_{n=0}^\infty \frac{(-1)^n \prod_{k=1}^n (p -2k +2) \prod_{j=1}^n (p+2j-1)}{(2n)!}x^{2n}$$
For $g_p$ where $p \in \{0,1,2,3\}$, we get the following solutions (see Problem 2 for more details):
\begin{align}
g_0 = 1\\
g_1 = -(1/9)x^10 - (1/7)x^8 - (1/5)x^6 - (1/3)x^4 - x^2 + 1\\
g_2 = -3x^2 + 1\\
g_3 = (2/7)x^10 + (3/7)x^8 + (4/5)x^6 + 3x^4 - 6x^2 + 1
\end{align}
If we set $a_0=0$ and $a_1=1$, we get $a_0 = a_2 = 0$. Thus, for $p \in \{0,1,2,3\}$ we get the following coefficients:
\begin{align}
a_3 = \frac{-(p+2)(p-1)}{3!}\\
a_5 = \frac{(p+4)(p+2)(p-3)(p-1)}{5!}
\end{align}
These initial conditions give us the independent solution $$h_p = a_1(x - \frac{(p+2)(p-1)}{3!}x^3 + \frac{(p+4)(p+2)(p-3)(p-1)}{5!}x^5 - ...) = a_1 \sum_{n=0}^\infty \frac{(-1)^n \prod_{k=1}^n (p -2k + 1) \prod_{j=1}^n (p+2j)}{(2n+1)!}x^{2n+1}$$
For $h_p$ where $p \in \{0,1,2,3\}$, we get the following solutions (see Problem 2 for more details):
\begin{align}
h_0 = (1/11)x^11 + (1/9)x^9 + (1/7)x^7 + (1/5)x^5 + (1/3)x^3 + x\\
h_1 = x\\
h_2 = -(2/33)x^11 - (5/63)x^9 - (4/35)x^7 - (1/5)x^5 - (2/3)x^3 + x\\
h_3 = -(5/3)x^3 + x
\end{align}

Chebyshev's equation (3) $$(1−x^2)y′′−xy′+p^2y=0$$
Suppose $y = \sum_{n  = 0}^\infty a_n x^n$ solves (3) Then we have $y' = \sum_{n  = 1}^\infty na_n x^{n-1}$ and $y'' = \sum_{n  = 2}^\infty (n-1)na_n x^{n-2}$.  Plugging these into (1), we get:
\begin{align}
\sum_{n  = 2}^\infty (n-1)na_n x^{n-2} - x^2\sum_{n  = 2}^\infty (n-1)na_n x^{n-2} - x\sum_{n  = 1}^\infty na_n x^{n-1} + p^2 \sum_{n  = 0}^\infty a_n x^n = 0  \\
\sum_{n  = 0}^\infty (n+1)(n+2)a_{n+2} x^{n} - \sum_{n  = 0}^\infty (n-1)na_n x^{n} - \sum_{n  = 1}^\infty na_n x^{n} + \sum_{n  = 0}^\infty p^2a_n x^n = 0 \\
\sum_{n  = 0}^\infty ((n+1)(n+2)a_{n+2} + (n-1)na_n - na_n + p^2a_n)x^n = 0 \\
(n+1)(n+2)a_{n+2} + (n-1)na_n - na_n + p^2a_n = 0
\end{align}
This gives the recurrence relation: $$a_{n+2} = \frac{(n-p)(n+p)}{(n+1)(n+2)}a_n$$
If we set $a_0=1$ and $a_1=0$, we get $a_1 = a_3 = a_5 = 0$. Thus we get the following coefficients:
\begin{align}
a_2 = \frac{p^2}{2!} \\
a_4 = \frac{(2^2 - p^2)p^2}{4!}\\
\end{align}
These initial conditions give us the independent solution $$g_p = a_0(1 + \frac{p^2}{2!}x^2 + \frac{(2^2 - p^2)p^2}{4!}x^4 + ...) = a_0 \sum_{n=0}^\infty \frac{\prod_{k=0}^n (p^2 -(2k)^2){(2n)!}x^{2n}$$

If we set $a_0=0$ and $a_1=1$, we get $a_0 = a_2 = 0$. Thus, for $p \in \{0,1,2,3\}$ we get the following coefficients:
\begin{align}
a_3 = \frac{1^2 - p^2}{3!}\\
a_5 = \frac{(3^2 - p^2)(1^2 - p^2)}{5!}
\end{align}
These initial conditions give us the independent solution $$h_p = a_1(x + \frac{1^2 - p^2}{3!}x^3 + \frac{(3^2 - p^2)(1^2 - p^2)}{5!}x^5 + ...) = a_1 \sum_{n=0}^\infty \frac{\prod_{k=0}^n (p^2 -(2k+1)^2){(2n+1)!}x^{2n+1}$$
The independent polynomial solutions are computed in the code below:


```python
# y_1 is even terms, we get this when a_0 = 1
# y_2 is odd terms, we get this when a_1 = 1
# y = a_0*y_1 + a_1*y_2

import math

x = var('x')


for p in range(0,4):

    g=1

    for n in range(1,3):

        prod1 = 1

        for k in range(0,n+1):

            prod1 = prod1*((2*k)^2 - p^2)

        g=g+((prod1/math.factorial(2*n))*x^(2*n))

    print("At p =",p,",\n","y_1 gives", g) 
    
    h=x

    for n in range(1,3):

        prod2 = 1

        for k in range(0,n+1):

            prod2 = prod2*((2*k + 1)^2 - p^2)

        h=h+((prod2/math.factorial(2*n+1))*x^(2*n+1))

    print(" y_2 gives", h) 

    print()
```

    At p = 0 ,
     y_1 gives 1
     y_2 gives 15/8*x^5 + 3/2*x^3 + x
    
    At p = 1 ,
     y_1 gives -15/8*x^4 - 3/2*x^2 + 1
     y_2 gives x
    
    At p = 2 ,
     y_1 gives 1
     y_2 gives -21/8*x^5 - 5/2*x^3 + x
    
    At p = 3 ,
     y_1 gives 105/8*x^4 + 45/2*x^2 + 1
     y_2 gives x
    


This doesn't look quite right, there must be a bug in my code.  It should be $y_1 = 1$ for $p=0$, $y_2 = x$ for $p = 1$, $y_1 = 2x^2-1$ for $p=2$, and $y_2 = 4x^3 -3x$ for $p=3$.

<h3>
    Problem 2
</h3>
Forp=0,1,2,3,4,5,6, find polynomial solutions for equation (2) with y(1)=1. Argue that, for each p there is only one such solution.  These polynomials are called the Legendre polynomials


```python
# y_1 is even terms, we get this when a_0 = 1
# y_2 is odd terms, we get this when a_1 = 1
# y = a_0*y_1 + a_1*y_2

import math

x = var('x')


for p in range(0,7):

    g=1

    for n in range(1,6):

        prod1 = 1
        
        prod3 = 1

        for k in range(1,n+1):

            prod1 = prod1*(p - 2*k + 2)
            
        for j in range(1,n+1):
            
            prod3 = prod3*(p+2*j-1)

        g=g+((-1)^n*(prod1*prod3/math.factorial(2*n))*x^(2*n))

    print("At p =",p,",\n","y_1 gives", g) 
    
    
    h=x

    for n in range(1,6):

        prod2 = 1
        
        prod4 = 1

        for k in range(1,n+1):

            prod2 = prod2*(p - 2*k + 1)
            
        for j in range(1,n+1):
            
            prod4 = prod4*(p+2*j)

        h=h+((-1)^n*(prod2*prod4/math.factorial(2*n+1))*x^(2*n+1))

    print(" y_2 gives", h) 

    print()
```

    At p = 0 ,
     y_1 gives 1
     y_2 gives 1/11*x^11 + 1/9*x^9 + 1/7*x^7 + 1/5*x^5 + 1/3*x^3 + x
    
    At p = 1 ,
     y_1 gives -1/9*x^10 - 1/7*x^8 - 1/5*x^6 - 1/3*x^4 - x^2 + 1
     y_2 gives x
    
    At p = 2 ,
     y_1 gives -3*x^2 + 1
     y_2 gives -2/33*x^11 - 5/63*x^9 - 4/35*x^7 - 1/5*x^5 - 2/3*x^3 + x
    
    At p = 3 ,
     y_1 gives 2/7*x^10 + 3/7*x^8 + 4/5*x^6 + 3*x^4 - 6*x^2 + 1
     y_2 gives -5/3*x^3 + x
    
    At p = 4 ,
     y_1 gives 35/3*x^4 - 10*x^2 + 1
     y_2 gives 1/11*x^11 + 1/7*x^9 + 2/7*x^7 + 6/5*x^5 - 3*x^3 + x
    
    At p = 5 ,
     y_1 gives -x^10 - 15/7*x^8 - 10*x^6 + 30*x^4 - 15*x^2 + 1
     y_2 gives 21/5*x^5 - 14/3*x^3 + x
    
    At p = 6 ,
     y_1 gives -231/5*x^6 + 63*x^4 - 21*x^2 + 1
     y_2 gives -8/33*x^11 - 5/9*x^9 - 20/7*x^7 + 10*x^5 - 20/3*x^3 + x
    


We can argue that the Legendre Equation has odd solutions for odd values of p and even solutions for even values of p, because we can see above that the degrees of the polynomials don't match up for p values and n that are not both even or both odd.  Thus, there is only one polynomial solution $L_p$ for each value of $p$: 
\begin{align}
L_0 = 1 \\
L_1 = x \\
L_2 = -3x^2 + 1 \\
L_3 = -\frac{5}{3}x^3 + x
\end{align}

<h3>
    Problem 3
</h3>
Plot the first seven Legendre polynomials on a single plot and conjecture aproperty for these polynomials that depends on whetherpis even or odd.


```python
# plot of the first 7 Legendre polynomials
var('x')
g = Graphics()
g += plot(1, (x,-1,1), ymin=-1,ymax=1,color='blue')
g += plot(x, (x, -1, 1), ymin=-1,ymax=1,color='red')
g += plot((.5*(3*x^2-1)), (x, -1, 1),ymin=-1,ymax=1, color='green')
g += plot((.5*(5*x^3-3*x)), (x,-1,1),ymin=-1,ymax=1, color='yellow')
g += plot(((1/8)*(35*x^3-30*x^2+3)), (x,-1,1), ymin=-1,ymax=1, color='teal')
g += plot(((1/8)*(63*x^5-70*x^3+15*x)), (x,-1,1), ymin=-1,ymax=1, color='pink')
g += plot(((1/16)*(231*x^6-315*x^4+105*x^2-5)), (x,-1,1), ymin=-1,ymax=1, color='purple')
g.show()
```

Based on this graph, it appears that the functions for odd p have odd order graphs and functions with even p have even order graphs.

<h3>
    Problem 4 
</h3>
Consider Hermite’s equation (1).  For a positive integer p, how should you choose values in $\{0,1\}$ for $a_0$ and $a_1$ so that the cooresponding power series solution $h_p(x)$ is actually a polynomial? Compute $h_p(x)$ for $p \in \{0,1,2,3,4,5,6\}$.


```python
# y_1 is even terms, we get this when a_0 = 1
# y_2 is odd terms, we get this when a_1 = 1
# y = a_0*y_1 + a_1*y_2

import math

x = var('x')


for p in range(0,7):

    g=1

    for n in range(1,7):

        prod1 = 1

        for k in range(1,n+1):

            prod1 = prod1*(2*k -2 - p)

        g=g+((2^n*prod1/math.factorial(2*n))*x^(2*n))

    print("At p =",p,",\n","y_1 gives", g) 
    
    
    h=x

    for n in range(1,7):

        prod2 = 1

        for k in range(1,n+1):

            prod2 = prod2*(2*k - 1 -p)

        h=h+((2^n*prod2/math.factorial(2*n + 1))*x^(2*n + 1))

    print(" y_2 gives", h) 

    print()
```

    At p = 0 ,
     y_1 gives 1
     y_2 gives 1/9360*x^13 + 1/1320*x^11 + 1/216*x^9 + 1/42*x^7 + 1/10*x^5 + 1/3*x^3 + x
    
    At p = 1 ,
     y_1 gives -1/7920*x^12 - 1/1080*x^10 - 1/168*x^8 - 1/30*x^6 - 1/6*x^4 - x^2 + 1
     y_2 gives x
    
    At p = 2 ,
     y_1 gives -2*x^2 + 1
     y_2 gives -1/102960*x^13 - 1/11880*x^11 - 1/1512*x^9 - 1/210*x^7 - 1/30*x^5 - 1/3*x^3 + x
    
    At p = 3 ,
     y_1 gives 1/23760*x^12 + 1/2520*x^10 + 1/280*x^8 + 1/30*x^6 + 1/2*x^4 - 3*x^2 + 1
     y_2 gives -2/3*x^3 + x
    
    At p = 4 ,
     y_1 gives 4/3*x^4 - 4*x^2 + 1
     y_2 gives 1/308880*x^13 + 1/27720*x^11 + 1/2520*x^9 + 1/210*x^7 + 1/10*x^5 - x^3 + x
    
    At p = 5 ,
     y_1 gives -1/33264*x^12 - 1/2520*x^10 - 1/168*x^8 - 1/6*x^6 + 5/2*x^4 - 5*x^2 + 1
     y_2 gives 4/15*x^5 - 4/3*x^3 + x
    
    At p = 6 ,
     y_1 gives -8/15*x^6 + 4*x^4 - 6*x^2 + 1
     y_2 gives -1/432432*x^13 - 1/27720*x^11 - 1/1512*x^9 - 1/42*x^7 + 1/2*x^5 - 5/3*x^3 + x
    


As you can see above, certain values for $a_0$ and $a_1$ give simpler polynomials.  For $p = 1$ I would choose $a_0 = 1, a_1 = 0$.  

<h3>
    Problem 5
</h3>
Let $G(x,t)=exp(2xt−t^2)$ and, for a positive integer $j$, define: $H_j(x)=\frac{_∂jG}{∂t^j}$
Show that $\frac{\partial^2G}{\partial x^2}−2x\frac{\partial G}{\partial x}+2(2xt−2t^2)G=0$.
The polynomials $H_n$ are called Hermite polynomials.  Make a conjecture for a formula that gives the scalars relating $h_n$ and $H_n$.
$\partial$


```python
# show that G'' - 2xG' + 2(2xt - 2t^2)G = 0
x,t = var('x,t')

G = e^(2*x*t - t^2)
dG = G.diff(x)
d2G = G.diff(x,2)

de = d2G - 2*x*dG + 2*(2*x*t - 2*t^2)*G
de.full_simplify()
```




    0




```python
# plug in some values of n to make a conjecture for scalars relating hn and Hn
Hlist = []
for n in range(0,7):
    new = G.diff(t,n)
    Hlist.append(new)
    print("H_",n,"=",new)
    
show(Hlist)
```

    H_ 0 = e^(-t^2 + 2*t*x)
    H_ 1 = -2*(t - x)*e^(-t^2 + 2*t*x)
    H_ 2 = 4*(t - x)^2*e^(-t^2 + 2*t*x) - 2*e^(-t^2 + 2*t*x)
    H_ 3 = -8*(t - x)^3*e^(-t^2 + 2*t*x) + 12*(t - x)*e^(-t^2 + 2*t*x)
    H_ 4 = 16*(t - x)^4*e^(-t^2 + 2*t*x) - 48*(t - x)^2*e^(-t^2 + 2*t*x) + 12*e^(-t^2 + 2*t*x)
    H_ 5 = -32*(t - x)^5*e^(-t^2 + 2*t*x) + 160*(t - x)^3*e^(-t^2 + 2*t*x) - 120*(t - x)*e^(-t^2 + 2*t*x)
    H_ 6 = 64*(t - x)^6*e^(-t^2 + 2*t*x) - 480*(t - x)^4*e^(-t^2 + 2*t*x) + 720*(t - x)^2*e^(-t^2 + 2*t*x) - 120*e^(-t^2 + 2*t*x)





<html></html>




```python
t = 0
H0 = e^(-t^2 + 2*t*x)
H1 = -2*(t - x)*e^(-t^2 + 2*t*x)
H2 = 4*(t - x)^2*e^(-t^2 + 2*t*x) - 2*e^(-t^2 + 2*t*x)
H3 = -8*(t - x)^3*e^(-t^2 + 2*t*x) + 12*(t - x)*e^(-t^2 + 2*t*x)
H4 = 16*(t - x)^4*e^(-t^2 + 2*t*x) - 48*(t - x)^2*e^(-t^2 + 2*t*x) + 12*e^(-t^2 + 2*t*x)
H5 = -32*(t - x)^5*e^(-t^2 + 2*t*x) + 160*(t - x)^3*e^(-t^2 + 2*t*x) - 120*(t - x)*e^(-t^2 + 2*t*x)
H6 = 64*(t - x)^6*e^(-t^2 + 2*t*x) - 480*(t - x)^4*e^(-t^2 + 2*t*x) + 720*(t - x)^2*e^(-t^2 + 2*t*x) - 120*e^(-t^2 + 2*t*x)

show("H_0 = G =",H0)
show("H_1 = dG =",H1)
show("H_2 = d2G =",H2)
show("H_3 = d3G =",H3)
show("H_4 = d4G =",H4)
show("H_5 = d5G =",H5)
show("H_6 = d6G =",H6)
```




<html></html>






<html></html>






<html></html>






<html></html>






<html></html>






<html></html>






<html></html>



From problem 4 we have $$h_0 = 1$$ $$h_1 = x$$ $$h_2 = -2x^2 + 1$$ $$h_3 = -\frac{2}{3}x^3 + x$$ $$\frac{4}{3}*x^4 - 4*x^2 + 1$$  and so on.  Thus, the equations for $h_n$ and $H_n$ both have order $n$ but are off by a coefficient.  I was unable to find a general pattern for to correct for the differences in the coefficients.

<h3>
    Problem 6
</h3>
Plot $h_j$ for $j=0,...,4$ on a single plot. Do the same for $H_j$ for $j=0,...,4$.


```python
var('x')
h = Graphics()
h += plot(1, (x,-2,2), ymin=-2,ymax=2,color='green')
h += plot(x, (x,-2,2), ymin=-2,ymax=2,color='red')
h += plot((-2*x^2 + 1), (x,-2,2), ymin=-2,ymax=2, color='blue')
h += plot((-2/3*x^3 + x), (x,-2,2), ymin=-2,ymax=2, color='teal')
h += plot((4/3*x^4 - 4*x^2 + 1), (x,-2,2), ymin=-2,ymax=2, color='magenta')
h.show()
```


```python
var('x')
H = Graphics()
H += plot(1, (x,-2,2), ymin=-2,ymax=2,color='green')
H += plot(2*x, (x, -2, 2), ymin=-2,ymax=2,color='red')
H += plot((4*x^2 - 2), (x, -2, 2),ymin=-2,ymax=2, color='blue')
H += plot((8*x^3 - 12*x), (x,-2,2),ymin=-2,ymax=2, color='teal')
H += plot((16*x^4-48*x^2+12), (x,-2,2), ymin=-2,ymax=2, color='magenta')
H.show()
```

<h3>
    Problem 7
</h3>
Show that $F(\alpha,\beta,\gamma,x)$ converges on $(−1,1)$ and is a solution for (5) correspondingto the indical root zero.

Let
\begin{align}
 P_0(x) & = -\alpha \beta, \\
 P_1(x) & = \gamma - (1+\alpha +\beta )x, \\
 P_2(x) & = x(1-x)
\end{align}
Then $$P_2(0) = P_2 (1)=0$$
We can see that x = 0 and x = 1 are singular points.  Let's use the ratio test to check if x = 0 is regular:
\begin{align}
\lim_{x \to a} \frac{(x - a) P_1(x)}{P_2(x)} &=\lim_{x \to 0} \frac{(x - 0)(\gamma - (1 + \alpha + \beta)x)}{x(1 - x)}=\lim_{x \to 0} \frac{x(\gamma - (1 + \alpha + \beta)x)}{x(1 - x)}= \gamma  \\ 
 \lim_{x \to a} \frac{(x - a)^2 P_0(x)}{P_2(x)} &= \lim_{x \to 0} \frac{(x - 0)^2(-\alpha \beta)}{x(1 - x)} = \lim_{x \to 0} \frac{x^2 (-\alpha \beta)}{x(1 - x)} = 0
\end{align}
Thus, x = 0 is a regular singular point and we can assume a solution of the form $$y = \sum_{n=0}^\infty a_n x^{n+m}$$  

Then, $$y' = \sum_{n=0}^\infty a_n(n+m)x^{n+m-1}$$ and $$y'' = \sum_{n=0}^\infty a_n(n+m)(n+m-1)x^{n+m-2}$$ Thus, after subbing these sums into the hypergeometric equation and simplifying we get the indicial equation:   $$a_{0} (m(m - 1) + \gamma m) = 0$$

Since $a_0 \not = 0$, we have the indicial roots $m = 0$ and $m = 1-\gamma$

The other terms give us the recurrence relation: $$a_n = \frac{(n + m + \alpha - 1)(n + m + \beta - 1)}{(n + m)(n + m + \gamma - 1)} a_{n - 1}, \text{ for } n \geq 1$$

After solving for a few terms, we can write the general formula: $$a_n =\frac{(m + \alpha)_n (m + \beta)_n}{(m + 1)_n(m + \gamma)_n} a_0, \text{ for } n \geq 0$$

This means our solution is of the form: $$y = a_0 \sum_{n = 0}^\infty \frac{(m + \alpha)_n (m + \beta)_n}{(m + 1)_n (m + \gamma)_n} x^{n + m}$$

Our formula for $a_n$ converges on(-1,1) by the ration test, the formula for $y$ we computed is a valid solution to the hypergeometric equation.

<h3>
    Problem 8
</h3>
Suppose that $\gamma$ is not an integer and express a solution to (5) corresponding to the indical root $1-\gamma$ in terms ofF. Give the general solution to (5) in terms of F.

If we base the particular solution $y_1$ on the indicial root $m = 0$ and $y_2$ on the indicial root $m = 1 - \gamma$, then we get:
\begin{align}
 y_1 &= a_0 \sum_{r = 0}^\infty \frac{(\alpha)_r (\beta)_r}{(1)_r (\gamma)_r} x^r = a_0 \cdot {{}_2 F_1}(\alpha, \beta; \gamma; x) \\ 
 y_2 &= a_0 \sum_{r = 0}^\infty \frac{(\alpha + 1 - \gamma)_r (\beta + 1 - \gamma)_r}{(1 - \gamma + 1)_r (1 - \gamma + \gamma)_r} x^{r + 1 - \gamma} \\
       &= a_0 x^{1 - \gamma} \sum_{r = 0}^\infty \frac{(\alpha + 1 - \gamma)_r (\beta + 1 - \gamma)_r}{(1)_r (2 - \gamma)_r} x^r \\ 
       &= a_0 x^{1 - \gamma} {{}_2 F_1}(\alpha - \gamma + 1, \beta - \gamma + 1; 2 - \gamma; x)
\end{align}

Thus our general solution is $y = A'y_1 + B'y_2$.  Also, if we let $A′ a_0 = a$ and $B′ a_0 = B$. Then $y = A_2F_1(\alpha ,\beta ;\gamma ;x) + Bx^{1-\gamma }_2F_1(\alpha -\gamma +1,\beta -\gamma +1;2-\gamma ;x)$

<h3>
    Problem 9
</h3>
Use the substitution $\frac{1}{2}(1−x)=t$ to show that a particular solution to Legendre’s equation (2) is given by $$F(p+1,−p,1,12(1−x))$$.


```python
t,x = var('t,x')
t = .5*(1-x)



for p in range(0,9):
    Alpha = p + 1
    Beta = -p
    Gamma = 1
    g=1

    for n in range(1,9):

        prod1 = 1
        
        prod3 = 1

        for k in range(1,n+1):

            prod1 = prod1*Alpha*Beta
            
        for j in range(1,n+1):
            
            prod3 = prod3*1*Gamma
            
        g=g+((prod1/prod3)*t^(n))

    print("At p =",p,",\n","y_1 gives", g.full_simplify()) 
    
    
    h=x

    for n in range(1,9):

        prod2 = 1
        
        prod4 = 1

        for k in range(1,n+1):

            prod2 = prod2*(Alpha + 1 - Gamma)*(Beta + 1 - Gamma)
            
        for j in range(1,n+1):
            
            prod4 = prod4*1*(2-Gamma)

        h=h+((prod2/prod4)*t^(n+1-Gamma))

    print(" y_2 gives", h.full_simplify()) 

    print()
```

    At p = 0 ,
     y_1 gives 1
     y_2 gives x
    
    At p = 1 ,
     y_1 gives 1.0*x^8 - 7.0*x^7 + 22.0*x^6 - 40.0*x^5 + 46.0*x^4 - 34.0*x^3 + 16.0*x^2 - 4.0*x + 1.0
     y_2 gives 1.0*x^8 - 7.0*x^7 + 22.0*x^6 - 40.0*x^5 + 46.0*x^4 - 34.0*x^3 + 16.0*x^2 - 3.0*x
    
    At p = 2 ,
     y_1 gives 6561.0*x^8 - 50301.0*x^7 + 169128.0*x^6 - 325620.0*x^5 + 392526.0*x^4 - 303318.0*x^3 + 146700.0*x^2 - 40596.0*x + 4921.0
     y_2 gives 6561.0*x^8 - 50301.0*x^7 + 169128.0*x^6 - 325620.0*x^5 + 392526.0*x^4 - 303318.0*x^3 + 146700.0*x^2 - 40595.0*x + 4920.0
    
    At p = 3 ,
     y_1 gives 1679616.0*x^8 - 13156992.0*x^7 + 45116352.0*x^6 - 88452000.0*x^5 + 108437616.0*x^4 - 85121064.0*x^3 + 41779836.0*x^2 - 11723034.0*x + 1439671.0
     y_2 gives 1679616.0*x^8 - 13156992.0*x^7 + 45116352.0*x^6 - 88452000.0*x^5 + 108437616.0*x^4 - 85121064.0*x^3 + 41779836.0*x^2 - 11723033.0*x + 1439670.0
    
    At p = 4 ,
     y_1 gives 100000000.0*x^8 - 790000000.0*x^7 + 2731000000.0*x^6 - 5395900000.0*x^5 + 6664510000.0*x^4 - 5269039000.0*x^3 + 2604057100.0*x^2 - 735537190.0*x + 90909091.0
     y_2 gives 100000000.0*x^8 - 790000000.0*x^7 + 2731000000.0*x^6 - 5395900000.0*x^5 + 6664510000.0*x^4 - 5269039000.0*x^3 + 2604057100.0*x^2 - 735537189.0*x + 90909090.0
    
    At p = 5 ,
     y_1 gives 2562890625.0*x^8 - 20332265625.0*x^7 + 70576312500.0*x^6 - 140001412500.0*x^5 + 173589378750.0*x^4 - 137762214750.0*x^3 + 68336450100.0*x^2 - 19371849060.0*x + 2402709961.0
     y_2 gives 2562890625.0*x^8 - 20332265625.0*x^7 + 70576312500.0*x^6 - 140001412500.0*x^5 + 173589378750.0*x^4 - 137762214750.0*x^3 + 68336450100.0*x^2 - 19371849059.0*x + 2402709960.0
    
    At p = 6 ,
     y_1 gives 37822859361.0*x^8 - 300781786347.0*x^7 + 1046518208442.0*x^6 - 2080767777480.0*x^5 + 2585828322126.0*x^4 - 2056717275354.0*x^3 + 1022463993096.0*x^2 - 290470182324.0*x + 36103638481.0
     y_2 gives 37822859361.0*x^8 - 300781786347.0*x^7 + 1046518208442.0*x^6 - 2080767777480.0*x^5 + 2585828322126.0*x^4 - 2056717275354.0*x^3 + 1022463993096.0*x^2 - 290470182323.0*x + 36103638480.0
    
    At p = 7 ,
     y_1 gives 377801998336.0*x^8 - 3008923058176.0*x^7 + 10484487344128.0*x^6 - 20876434539520.0*x^5 + 25981030302976.0*x^4 - 20694127547968.0*x^3 + 10302164328400.0*x^2 - 2930773171396.0*x + 364774343221.0
     y_2 gives 377801998336.0*x^8 - 3008923058176.0*x^7 + 10484487344128.0*x^6 - 20876434539520.0*x^5 + 25981030302976.0*x^4 - 20694127547968.0*x^3 + 10302164328400.0*x^2 - 2930773171395.0*x + 364774343220.0
    
    At p = 8 ,
     y_1 gives 2821109907456.0*x^8 - 22490515095552.0*x^7 + 78444705042432.0*x^6 - 156349507599360.0*x^5 + 194767298862336.0*x^4 - 155282346730944.0*x^3 + 77377486975056.0*x^2 - 22033095055164.0*x + 2744863693741.0
     y_2 gives 2821109907456.0*x^8 - 22490515095552.0*x^7 + 78444705042432.0*x^6 - 156349507599360.0*x^5 + 194767298862336.0*x^4 - 155282346730944.0*x^3 + 77377486975056.0*x^2 - 22033095055163.0*x + 2744863693740.0
    


These polynomials resemble the Legendre polynomials.  We already proved that these are solutions to the hypergeometric equation.

<h3>
    Problem 10: Curiosity Element
</h3>
I chose to investigate the differential equations associated with the SIR model infectious disease spread from my Systems Ecology Modelling independent study last spring.  We can apply the model to the COVID 19 pandemic.

More info about the model I used can be found here: https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model, https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-contact-number

To model the global spread of COVID-19 I will define 3 differential equations.  To start define  $t$ measured in days, and 3 dependent variables: $S = S(t)$ is the number of susceptible individuals, $I = I(t)$ is the number of infected individuals, and $R = R(t)$ is the number of recovered individuals.  We can define 3 more dependent variables: $s(t) = S(t)/N$ is the proportion of susceptible individuals in a population of $N$ individuals, $i(t) = I(t)/N$ is the proportion of the total population that is infected, and $r(t) = R(t)/N$ is the proportion of the total population that is recovered.  This model assumes that people cannot get Sars-CoV-2 more than once, which is not neccessarily true.  We can now define three ordinary differential equations:
\begin{align}
\frac{ds}{dt} = -bs(t)i(t) \\
\frac{dr}{dt} = ki(t) \\
\frac{di}{dt} = bs(t)i(t) - ki(t)
\end{align}
Where $b$ is a fixed number of contacts per day with infected individuals that applies to every individual, and $k$ is a fixed fraction of the infected group that recovers per day.

We know that a typical COVID patient is infectious for about 14 days, so we can set $k = 1/14$.  Let $b = 1/2$, meaning each person comes in contact with an infectious COVID patient every 2 days.  Now we can try to solve and plot the equations with sage.


```python
# here t was changed to x to check for code errors
from sage.symbolic.integration.integral import indefinite_integral

x = var('x')
s = function('s')(x)
i = function('i')(x)
r = function('r')(x)

sus = desolve(diff(s,x) + (1/3)*i*s,s)
imm = desolve(diff(i,x)-(1/3)*s*i + (1/14)*i,i)
rec = desolve(diff(r,x)-(1/14)*i,r)

print("s = ")
sus.show()
print("i = ")
imm.show()
print("r = ")
rec.show()

rec = (1/14)*indefinite_integral(imm,x)
print("r =")
rec.show()
```

    s = 





<html></html>



    i = 





<html></html>



    r = 





<html></html>



    r =





<html></html>



I am unable to graph the equations because they all depend on each other.  In order to graph I will need to do more research on solving systems of differential equations. This technique may involve the Wronskian and variation of parameters.  I don't have time to try this now, so maybe another time. 


```python

```

af58cdde28ad97c873c1bbc958f77dbb82b0f8d8
ipynb
Jupyter Notebook
LyonDelSordo ODEfinalProject.ipynb
esther-lyondelsordo/Ordinary-Differential-Equations
9114489d638075a497ac52e1cb10714f871d7fa1
MIT
LyonDelSordo ODEfinalProject.ipynb
esther-lyondelsordo/Ordinary-Differential-Equations
9114489d638075a497ac52e1cb10714f871d7fa1
MIT
LyonDelSordo ODEfinalProject.ipynb
esther-lyondelsordo/Ordinary-Differential-Equations
9114489d638075a497ac52e1cb10714f871d7fa1
MIT
Qwen/Qwen-72B
1. YES
2. YES