共计 4712 个字符，预计需要花费 12 分钟才能阅读完成。
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.
- 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. - 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 - 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 - 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. - 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.
http://www.daixie0.com/conten…
正文完