ass2留学生代做、代写Python程序设计、Python语言代做、programs留学生代写

日期：2019-05-30 11:03

ass2

May 29, 2019

1 ASSIGNMENT 2: SUMS & DIFFERENCES

DUE: 11PM 3 June, 2019

1.0.1 INSTRUCTIONS:

0) Please cite any material (code or information) you use from the Internet. You can do this in

comments in your code, or separate Markdown cells.

1) Please feel free to work in groups. It’s more fun and more effective. However, if you do

work in groups, you must state the people you worked with here:

TEAM MEMBERS:

In [ ]: ### LIST TEAM HERE ###

EVERY PERSON IN THE GROUP MUST SUBMIT THEIR OWN ASSIGNMENT

EVERY PIECE OF CODE MUST SAY WHO IN THE GROUP WROTE IT. OR WHERE IT

CAME FROM IF YOU GET IT FROM SOMEWHERE ELSE. YOU NEED TO ALSO PUT THIS

INFORMATION IN THE COMMENTS.

2) Make sure to save your work. I recommend working out of a new notebook that you

create with different name than this file. This will make it so you’ll never lose any work if I

make an update at some point.

3) I will make sure your code runs without errors and gives the correct results in a reasonable

amount of time. I do not care about the fine details of how you get it to work.

4) The programs in Problems 1a-b, and Problem 2 do not need each other to work.

5) All the required material for this Assignment has been covered in the Lectures and Tutorials.

You can find examples in the online notebooks.

6) You DO NOT need to output any plots for you assignment. But you will probably want

to make your own plots to see what’s going on.

7) READ EVERYTHING CAREFULLY.

8) GOOD LUCK!

In [ ]: import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

1

2 Problem 1:

2.1 1a: Guass-Legendre quadrature:

Consider the three-term recursion relation

Pn1(x), with P0(x) = 1, P1(x) = x

This generates a sequence called Legendre polynomials. For example,

Now consider a second three-term recursion relation.

Qn+1(x) = (2 n + 1) Pn(x) + Qn1(x), with Q0(x) = 0, Q1(x) = 1

It possible to prove that

This means that Pn(x) is an nth-order polynomial, and Qn(x) is an (n ? 1)th-order polynomial.

2.1.1 TASK 1:

Make a function that takes a numpy array of x values (of any length) and returns all the polynomals

Pn(x) and Qn(x) up to order n=N. Hint: Fill in the given zero arrays with the higher-order P[n] and

Q[n].

In [ ]: ### DO NOT ERASE THIS CELL ###

def Legendre(N,x):

P = np.zeros((N+1,len(x)))

Q = np.zeros((N+1,len(x)))

P[0], P[1] = 1 + 0*x, x

### PUT CODE HERE ###

return P, Q

In [ ]: ## TEST CODE HERE ##

2

We are interested in the N roots of the PN. That is,

PN( x i ) = 0

The roots are not known exactly. But they are close to

If we know the roots, then we can define a set of set of N weights

Write a function called quadrature that computes the roots x, and weights, w, and returns the

results as length-N numpy arrays.

Hint: you’ll have to use Newton’s method for this. This is why you need to compute both P[N] and

Q[N]. Use the above formula as an initial guess. You can monitor the error and/or number of iterations how

ever you want. Just make sure the results are as converged as possible.

In [ ]: ### DO NOT ERASE THIS CELL ###

def quadrature(N,imax=10):

i = np.arange(N)

x = np.cos(np.pi*(4*i+3)/(4*N+2))

## PUT CODE HERE ##

## RETURN THE FOLLOWING:

return x, w

In [ ]: ## TEST CODE HERE ###

2.1.3 TASK 3:

It’s a fact that

Find (guess) the (simple) analytical form of μ(x). Hint: You should be able to deduce the shape of

μ(x) from a plots of w versus x.

Make a function that returns μ(x). The function should return an array of values if x s an array.

In [ ]: ### DO NOT ERASE THIS CELL ###

def mu(x):

3

## PUT CODE HERE ##

return

2.1.4 Guass-Legendre Integration:

The whole point of computing x and w is the following:

This integration method will converge much more rapidly than the standard Simpson’s etc. In

fact it’s exact if f(x) is any polynomial up to degree 2N + 1.

From now on, everything is easier if we always assume N is even.

2.1.5 TASK 4:

Make a function called integrate that accepts an array y, the weights w, and returns the integral.

The y array will be the result of computing y=f(x) for some function, f(x).

In [ ]: ### DO NOT ERASE THIS CELL ###

def integrate(y,w):

## PUT CODE HERE ##

In [ ]: ## TEST CODE HERE ###

2.1.6 TASK 5:

Compute the integral using your quadrature scheme.

Write a function that produces an array with the relative error in computing the above integral

for 4 ≤ n ≤ N, where n and N are even.

In [ ]: ### DO NOT ERASE THIS CELL ###

def relative_error(N):

### PUT CODE HERE ###

In [ ]: ### TEST CODE HERE ###

4

2.1.7 TASK 6:

Do the same thing as TASK 5, only using the trapezoidal rule.

In [ ]: ### DO NOT ERASE THIS CELL ###

def relative_trap_error(N):

### PUT CODE HERE ###

In [ ]: ### TEST CODE HERE ###

2.2 1b: Integrating when we don’t even really know the function.

Consider the function, y(t) given by the implicit relation

There is no closed-form solution for y(t). But clearly y = 0 is a solution when t = 0. And both y and t increase together. This means that this is a one-to-one relationship and we can solve

it numerically.

Now consider the integral

2.2.1 TASK 7:

The quadrature scheme assumes the integration goes from ?1 to +1. To compute I(z) you’ll need

to transform from x to t and rescale w to dt.

Make a function that does the transformation and returns a new grid and weights, t, and dt:

In [ ]: ### DO NOT ERASE THIS CELL ###

def scale(z,x,w):

## PUT CODE HERE ##

return t, dt

In [ ]: ### TEST CODE HERE ###

2.2.2 TASK 8:

Make a function that takes t values and uses Newton’s method to compute the solution to the

equation

y e y t = 0

5

Hint: The initial guess y = log(t + 1) works well enough.

In [ ]: ### DO NOT ERASE THIS CELL ###

import numpy as np

def y(t):

# Initial guess

y = np.log(t+1)

## PUT CODE HERE ###

return y

In [ ]: ### TEST CODE HERE ###

2.2.3 TASK 9:

Given an array of t values, make function that returns

f(t) = sin(t)

1 + t y(t)

Hint: you should be able to use y(t) just like any other function.

In [ ]: ### DO NOT ERASE THIS CELL ###

def f(t):

## PUT CODE HERE ##

In [ ]: ### TEST CODE HERE ###

2.2.4 TASK 10:

Make a function called I that takes z, and a number of grid points, N and computes the integral

outlined in the introduction.

In [ ]: ### DO NOT ERASE THIS CELL ###

def I(z,N):

## PUT CODE HERE ##

return

In [ ]: ### TEST CODE HERE ###

6

3 Problem 2: Solving nonlinear ODEs:

3.1 2a: Initial-Value Problem:

Suppose we have a 2nd-order nonlinear ODE:

We can think of this as an initial value problem, and set initial conditions,

x(0) = a, y(0) = b.

We can then integrate this up until some later time, T, and ask about.

This kind of problem requires a time-stepping method.

3.1.1 TASK 1:

Solve the following initial-value problem using a 4th-order Runge-Kutta method.

Make a function that solves the above initial-value problem. The function should take a numpy

array of time values, initial conditions x0, y0, and the parameter α. The function should return the

values of x(t) and y(t) at the values of the time array.

Hints:

You can use np.diff to find time-step sizes from your time array.

You can define functions inside of functions.

You should not need to use any fancy classes. It’s simpler to copy code from lecture and

tutorial examples.

My solution took 11 lines. I can imagine many reasonable solutions taking up to ~30 lines.

But much more than this probably means you are making things much more complicated

than you need to.

I am including an example plot of what the solution should look like for one set of parameters.

I made the plot using:

t = np.arange(0,16,1/4096)

x, y = IVP(y,0,1,1)

The above input took about three seconds to compute the result.

7

It’s a good sign if you can reproduce the plot. But that is still no guarantee everything is

correct. Make sure you test things carefully.

I am going to test you code against these parameters, and others that I will choose at the

time.

In [ ]: import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

In [ ]: ### DO NOT ERASE THIS CELL ###

def IVP(t,x0,y0,alpha):

### PUT CODE HERE ###

return x, y

In [ ]: ## TEST CODE HERE ###

## This should work ##

t = np.arange(0,16,1/4096)

x, y = IVP(t,0,1,10)

3.1.2 TASK 2:

The for the above IVP, the following quantity should stay constant in time:

Make a function that returns the maximum change in Z for a solution to the above IVP for a

given α.

In [ ]: ### DO NOT ERASE THIS CELL ###

def max_change_Z(x,y,alpha):

### PUT CODE HERE ###

In [ ]: ## TEST CODE HERE ###

3.1.3 TASK 3:

Write a function called spectrum that takes arrays with t and x(t), and returns the range of frequencies,

along with the real fast Fourier transform (FFT) of x(t).

8

Hints:

You can use you function to make a log-log periodograms of each spectrum.

You should use numpy or scipy functions to give you the frequencies and amplitudes.

Using the right numpy functions, my code took 3 lines.

I made the example periodogram plotting the results of the following input

t =np.arange(0,128,1/4096)

x, y = IVP(t,0,1,10)

f, A = spectrum(t,x)

The above code took a few seconds to compute x,y, but everything else is fast.

In [ ]: ### DO NOT ERASE THIS CELL ###

def spectrum(t,x):

## PUT CODE HERE ##

return frequency, amplitude

In [ ]: ## TEST CODE HERE ###

## This should work ##

t = np.arange(0,128,1/4096)

x, y = IVP(t,0,1,10)

f, A = spectrum(t,x)

3.2 2b: Boundary-Value Problem:

We can have the same ODE as Problem 2a, but have a very different way of thinking about it.

3.2.1 TASK 4:

Write a function that solves the following non-linear 2-point boundary-value problem using Newton’s

method and 2nd-order finite differences.

α u(x) = 1 for 0 < x < 1, and u(0) = u(1) = 0

Your function should take a numpy array as an initial guess, and an argument for α, and speci-

fied number of iterations. It should return the updated guess and a Cauchy error estimate.

9

Hints:

Recall that if we want to solve the nonlinear problem

u′′ ( x ) = f( u(x) ) where u(0) = u(1) = 0.

Then we need to iteratively solve:

( uold(x) ) uold(x) where unew(0) = unew(1) = 0.

You have to do is setup and solve a sequence of linear 2-point boundary value problems.

For any given iteration, you can define the Cauchy error

error = max

| unew(x) uold(x)|

As before, you shouldn’t need to use any fancy classes.

My solution took 14 lines. I can imagine many reasonable solutions taking up to ~40 lines.

But much more than this probably means you are making things much more complicated

than you need to.

I am including an example plot of what the solution should look like for one set of parameters.

I made the plot using:

x = np.linspace(0,1,50)

u, error = BVP(16*x*(1-x),1,iterations=5)

The result produced the plot, and the error was 3.323955022098061e-14

The above input took about less than one second to compute the result.

It’s a good sign if you can reproduce the plot. But that is still no guarantee everything is

correct. Make sure you test things carefully.

I am going to test you code against these parameters, and others that I will choose at the

time.

In any case I test, I will put in a guess that I know will converge quickly.

In [ ]: ### DO NOT ERASE THIS CELL ###

def BVP(guess,alpha,iterations=5):

### PUT CODE HERE ###

10

# First make 2nd-order 2nd derivative matrix.

# Then do Newton iterations.

return u, error

In [ ]: ### TEST CODE HERE ###

# This should work #

x = np.linspace(0,1,50)

u, error = BVP(16*x*(1-x),1,iterations=5)

4 END ASSIGNMENT

In [ ]:

11