代写EN 560.601、MATLAB编程代做、代写MATLAB设计、代做column vector

日期：2019-04-24 10:42

EN 560.601 - Spring 2019

Homework #9

Due on 12PM April 24th, 2019

This problem set contains significant portion of coding. So please attach

your code and the result in your HW.

1. Using Fourier Integral to obtain the solution of heat equation ut = c

2uxx in (∞, +∞), satisfying the following initial condition u(x, 0) = f(x).

Please verify that u(x, t) satisfy the initial condition.

(a) f(x) = 1/(1 + x2)

(b) f(x) = sin(x)/x, Please plot the solution u(x, t) with c = 1.

2. Step size for Differentiation

Consider the polynomial function f(x) in the domain x ∈ [0, 0.3]

We can calculate the first derivative of this polynomial analytically.

(a) Discretize the domain [0, 0.3] first using step size 0.1, (xk=0:0.1:0.3).

Use the O(x2) accurate center difference scheme to find out

the values of the first derivative f0(xk) at these discretized points.

(At boundaries use O(x2) accurate forward- and backward-schemes.)

Store these values as a column vector. Compute the maximum

of the absolute error over all values of f0(xk).

(b) Repeat the same problem with different step sizes

. Find the maximum of the absolute errors over

all values of f0

(xk) in a 1 × 6 row vector. Please plot these errors

using a logarithmic scale for the y-axis. In Matlab, the command

semilogy might help. From this figure, you will find out there

is one step size that has the smallest absolute error. Find this

optimal step size xopt. What is the theoretical optimal step size?

Does the numerical result fit the theoretical result?

1

3. Determining the Global Error of an Numerical Integration

Scheme

Given an numerical integration scheme, e.g. the Trapezoidal Rule and

the Simpson’s Rule, it is generally difficult to derive the global truncation

error analytically. In the case of the Trapezoidal Rule, the global

truncation error is

where a and b are left and right boundaries of the integral and c ∈ [a, b].

We want to verify such an error term of an arbitrary integration scheme

numerically.

Suppose an engineer proposed a new integration method and claimed

it is very accurate. In the first three subintervals of the integration

domain [x0, x3],

So the integration formula for the whole domain [a, b] is

Z xN

f(x0) + 3f(x1) + 3f(x2) + 2f(x3) + 3f(x4)

+ 3f(x5) + 2f(x6) + . . . + 3f(xN?1) + f(xN )

(1)

In order to implement this method, N should have the form of N =

3k, k = 1, 2, . . ..

Your boss wants you to verify this method and find out how accurate

this method could be.

(a) First we want to determine the order of accuracy for this method.

We pick a test function f(x) = exp(7x) on the domain [0, 1.5].

The antiderivative of this function can be easily calculated, that

is F = exp(7x)/7. Pick N = 399 and construct equidistant points

2

for x. There should be N + 1 points in total. Implement this

numerical integral method on the domain [0, 1.5] in MATLAB

and compare with true value (Qtrue=F(1.5)-F(0)). Find the

absolute value of the error (err1=abs(Q1-Qtrue)).

Repeat the previous step for N = 798. You can see ?x is almost

decreased by half now. Find the absolute value of the new error

(err2=abs(Q2-Qtrue)). Let’s divide err1 by err2 and see the

factor the global error changes when ?x is decreased by half.

Based on this, you are able to determine the order of accuracy for

this method.

(b) Second we are interested in the derivative term of the global error.

For the Trapezoidal Rule, the derivative term is the second

derivative of f(x). This implies if f(x) is a linear function, say,

f(x) = ax + b, the global error is always 0. Because f

00(x) is always

0. We will use this idea to find out which derivative term

appears in the error.

Use N = 399 and construct equidistant points for x. But this time

change the test function to six polynomials with different degrees.

Here choose f(x) to be x

. Again, implement this

numerical integration method on the domain [0, 1.5] in MATLAB

for these six different f(x) and compare with true values. Find the

absolute values of error for each f(x) into a 1 × 6 row vector one

by one. Based on the error values, you can tell which derivative

this error term has. Some values are about 10?16, and you can

take these values as 0 since they are near the round-off error.

(c) You can almost determine the global error term of this method.

The only part left is the coefficient term. In fact, you can find out

the coefficient based on the result in the second part. Try it out.

4. Van Der Pol Oscillator

Consider the van der Pol oscillator

In order to integrate this equation with the standard solution techniques,

it has to be written as a system of first order differential equations.

(a) Solve the equation for t = [0 : 0.01 : 32] using ode45. The initial

conditions are y(t = 0) = 2 and dy(t = 0)/dt = 0. Plot the

solution and note the intermittent behavior.

3

(b) We want to investigate how many steps ode45 requires. Solve the

same problem, but now do not specify to ode45 where to evaluate

the solution (by only supplying the start and end times [0, 32]).

Also, solve the same problem with ode15s. Save the number of

time-steps for both of methods as a 2×1 vector with ode45 results

first and ode15s results second.

Notes: Plot both solutions with points instead of solid lines. What do

you see? ode15s solver is appropriate for stiff ODEs. This means that

the solver can adjust its step-size to be small when accuracy is needed in

regions of rapid changes and to be bigger in regions where the solution

is changing slowly without worrying about having to maintain stability.

Even though the cost per time step is generally more expensive than

a non-stiff solver, the stiff solvers usually are able to take bigger timesteps

on average, and as a result solve the problem quicker.

4