Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

MATH 157 final project

Views: 84
Kernel: SageMath 8.1

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

Numerical Solutions to ODE

Preview of Ordinary Differential Equations

An ordinary differential equation (frequently called an "ODE") is an equality involving a function and its derivatives. An ODE of order n is an equation of the form F(x,y,y′,...,y(n))=0F(x,y,y',...,y^{(n)})=0 where y is a function of x, yny^{n} is the nthn^{th} derivative of x.

Sovling Differential Equations with Sage

Example: Solve x′+x=5,x(0)=1x'+x = 5, x(0) = 1

t = var('t') x = function('x', t) de = lambda y: diff(y,t) + y - 5 x0 = 1 soln = desolve(de(x),[x,t],[0,x0]) show(soln)
(5 et−4)e(−t)\renewcommand{\Bold}[1]{\mathbf{#1}}{\left(5 \, e^{t} - 4\right)} e^{\left(-t\right)}

Example: Solve x′′−x′=0x''-x' = 0

t = var('t') de = lambda x: diff(x,t,t) - x x0 = 1 x1 = 0 desolve(de(x),[x,t])
_K2*e^(-t) + _K1*e^t

Example: Solve for y in terms of t in 12(ln(y−1)−ln(y+1))=t+C\frac{1}{2}(ln(y-1) - ln(y+1)) = t+C

C = var('C') y = var('y') solve([log(y - 1)/2 - (log(y + 1)/2) == t + C],y)
[log(y + 1) == -2*C - 2*t + log(y - 1)]

Note we did not get an expression of y directly; let's try a different input.

solve([log((y - 1)/(y + 1)) == 2*t + 2*C],y)
[y == -(e^(2*C + 2*t) + 1)/(e^(2*C + 2*t) - 1)]

Solve the same equation by adding an initial condition: y(1) = 3

solny=lambda t:(-e^(2*C + 2*t)-1)/(e^(2*C + 2*t)-1) solve([solny(1) == 3],C)
[C == -1/2*log(2) - 1]
C = -1/2*log(2) - 1 solny(t)
-(e^(2*t - log(2) - 2) + 1)/(e^(2*t - log(2) - 2) - 1)

Visualizing Solutions to ODE

The idea of direction field (or slope field) related to the first order ODE is similar to vector calculus.

At each point (x,y)(x,y), we plot a small vector that has slope f(x,y)f(x,y). A collection of all points (x,y)(x,y) where the vectors have slope m generate direction field.

Example: Newton's Equation (MATLAB assignment exercise 2.5)

Newton's law of cooling models the temperature change of an object at a certain temperature when placed in a surrounding environment of a different temperature. The law can be stated as follows: dy/dt=k(A−y)dy/dt = k(A - y) where k is a positive proportionality constant and A represents temperature of the environment.

Plot a direction field with A = 1, k = 2, and where the minimum value of t is zero.(We start range of t from 0 because we do not need to worry about negative times)

t,y=var('t,y') f(t,y) = 2-2*y v=plot_slope_field(f(t,y),(t,0,10),(y,-10,10),headaxislength=3, headlength=3) show(v)
Image in a Jupyter notebook

We can see that all solutions approach to 1 which is environmental temperature under our setting.

Euler's Method

What if the differential equation cannot be solved with a nice formula? We introduce Euler's Method here. Euler's Method is an algorithmic way of plotting an approximate solution to an initial value problem through the direction field.

A more detailed explanation to Euler's method could be found here.

Recall from Calculus that f(x,y)≃y(x+h)−y(x)hf(x,y)\simeq \frac{y(x+h)-y(x)}{h}, solving for y(x+h)y(x+h) we get: y(x+h)≃y(x)+h⋅f(x,y(x))y(x+h)\simeq y(x)+h\cdot f(x,y(x)) where h is the step size.

Note: the first order DE must be in the form of y′=f(x,y),y(a)=cy'=f(x,y), y(a)=c for this method to work.

Example: Use Euler’s method with h=1/2h = 1/2 to approximate y(1)y(1), where y′−y=4x−4,y(0)=1y'-y=4x-4, y(0)=1

x,y=PolynomialRing(QQ,2,"xy").gens() eulers_method(4*x+y-4,1,1,1/3,2) eulers_method(4*x+y-4,0,1,1/2,1) pts = eulers_method(4*x+y-4,0,1,1/2,1,method="none") list_pts = list_plot(pts) line_pts = line(pts) show(list_pts+line_pts)
x y h*f(x,y) 1 1 1/3 4/3 4/3 8/9 5/3 20/9 44/27 2 104/27 212/81 x y h*f(x,y) 0 1 -3/2 1/2 -1/2 -5/4 1 -7/4 -7/8
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:5: DeprecationWarning: use the option 'algorithm' instead of 'method' See http://trac.sagemath.org/6094 for details.
Image in a Jupyter notebook

Further references

In order to improve acccuracy in Euler's method, we usually choose small step sizes. However, it is hard to avoid roundoff errors. If the step size is too small, we will build up so many roundoff errors into approximation. To solve the problem, we introduce an improved Euler's method, which is also called Heun's Method. The "improved" tangent line approximation at (a,c)(a,c) is:

y(a+h)≃c+h⋅m+m′2=c+h⋅f(a,c)+f(a+h,c+h⋅f(a,c))2y(a+h)\simeq c+h\cdot\frac{m+m'}{2} = c+h\cdot\frac{f(a,c)+f(a+h,c+h\cdot f(a,c))}{2}

Solving improved Euler's method using sage.

Fourth-Order Runge-Kutta Method

With Euler's method, we are able to solve most cases. However, sometimes we encounter difficult cases that require more sophiscated methods. Thus, we introduce fourth-order Runge-Kutta method. This method involves computing four slopes and taking an average of them. A more detailed explanation could be found here.

Example: y′=−y+cos(x)y' = -y + cos(x)

x = var('x') y = function('y')(x) DE = diff(y,x) == -y + cos(x) desolve_rk4(DE,y,ics=[0,1],end_points=[0,n(2*pi)],step=0.1,output='slope_field')
Image in a Jupyter notebook

Exercises

1.Use Euler’s method to approximate x(1)x(1), where x′′−3x′+2x=1,x(0)=0,x′(0)=1x''-3x'+2x=1, x(0)=0, x'(0)=1

Solution: Let x1=x,x2=x′x_1 = x, x_2 = x', we trasform given equations to standard form: x1′=x2,x2′=1−2x1+3x2,x1(0)=0,x2(0)=1x_1' = x_2, x_2' = 1-2x_1+3x_2, x_1(0)=0, x_2(0)=1.

Take h=(1−0)3=13h = \frac{(1-0)}{3} = \frac{1}{3}

RR = RealField(sci_not=0, prec=4, rnd='RNDU') t, x, y = PolynomialRing(RR,3,"txy").gens() f = y g = 1-2*x+3*y L = eulers_method_2x2(f,g,0,0,1,1/3,1,method="none") eulers_method_2x2(f,g, 0, 0, 1, 1/3, 1) P1 = list_plot([[p[0],p[1]] for p in L]) P2 = line([[p[0],p[1]] for p in L]) show(P1+P2)
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/repl/ipython_kernel/__main__.py:5: DeprecationWarning: use the option 'algorithm' instead of 'method' See http://trac.sagemath.org/6094 for details.
t x h*f(t,x,y) y h*g(t,x,y) 0 0 0.35 1 1.4 1/3 0.35 0.88 2.5 2.8 2/3 1.3 2.0 5.5 6.5 1 3.3 4.5 12. 11.
Image in a Jupyter notebook

2.Using the following function definition to complete Runge-Kutta method and Euler method(Note these two methods use different ways to calculate increment)

def compare_methods(xstart,ystart,xfinish,f,nsteps = 10,tol = 10^(-5.0)): sol = [ystart] xvals = [xstart] h = N((xfinish-xstart)/nsteps) #your code goes here

Solution

def compare_methods(xstart,ystart,xfinish,f,nsteps = 10,tol = 10^(-5.0)): sol = [ystart] xvals = [xstart] h = N((xfinish-xstart)/nsteps) for step in range(nsteps): k1 = f(xvals[-1],sol[-1]) k2 = f(xvals[-1] + h/2,sol[-1] + k1*h/2) k3 = f(xvals[-1] + h/2,sol[-1] + k2*h/2) k4 = f(xvals[-1] + h,sol[-1] + k3*h) ek2 = f(xvals[-1]+h,sol[-1]+k1*h) #Runge-Kutta sol.append(sol[-1] + h*(k1+2*k2+2*k3+k4)/6) #Euler's method sol.append(sol[-1] + h*(k1+ek2)/2) xvals.append(xvals[-1] + h) return zip(xvals,sol)