1. Consider the following nonlinear system about xi, i = 1, 2, 3, 4:∫ x2
x1
sin(x3t
2 + x4)dt = −2.886389616512186e− 01,∫ x2
x1
cos(x3t
2 + x4)dt = 3.070929941495893e− 01,∫ x2
x1
tcos(x3t
2 + x4)dt = 6.104465468269571e− 01,
sin(x3x2 + x4) = −9.764160102906496e− 01.
(a) Implement a Matlab function for f(x) and a Matlab function for its Jacobian Jf (x).
In these Matlab function, use the composite Gauss-Legendre quadrature formula
with 100 subintervals and 2 local nodes to compute all the integrals involved.
Then, compute f(x(0)) and Jf (x
(0)) with
x(0) = [0.9; 2.5; 1.35; 1.6]
and present the numerical results in the following table
f1(x
(0))
f3(x
(0))
Jf (x
(0))(1, 1)
Jf (x
(0))(2, 3)
Jf (x
(0))(4, 1)
1
yehong
附注
(b) Solve f(x) = 0 by Newton’s method started from x(0) = [0.9; 2.5; 1.35; 1.6]. Use
tol = 10^(-10); nmax = 200;
Present your solution in the following table:
x1
x2
x3
x4
num of iterations used
2. The top profile of a noble beast is the graph of a piecewise function f(x) digitized into
three data sets as follows:
Piece 1:
i XXXXXXXXXX
xi XXXXXXXXXX17
f(xi XXXXXXXXXX XXXXXXXXXX
f ′(xi XXXXXXXXXX
Piece 2:
i XXXXXXXXXX
xi XXXXXXXXXX27.7
f(xi XXXXXXXXXX XXXXXXXXXX
f ′(xi XXXXXXXXXX
Piece 3:
i XXXXXXXXXX
xi XXXXXXXXXX
f(xi XXXXXXXXXX
f ′(xi XXXXXXXXXX
Instruction: Use the composite Gauss-Legendre formula with 50 subintervals and 3
local nodes to compute integrals in this problem.
(a) Make a plot for the top profile of the noble beast by the clamped cubic spline.
Scale the plot by
2
axis([0,30,0,30])
In addition, use cubic spline interpolants to find values of f(x) at x values specified
y the following table and fill in the table the numerical results:
x f(x) by the c-CS f(x) by the n-CS
10.5
20.5
28.5
where “c-CS” means the cubic spline with the clamped boundary condition and
“n-CS” means the cubic spline with the natural boundary condition.
(b) Find the length L for the 2nd piece of the graph of f(x) in the top profile of the
noble beast.
(c) In this part, we cut the 2nd piece of the graph of f(x) into three equal pieces.
To do this, find x1, x2 ∈ [17, 27.7] such that the lengthes of the 2nd piece of the
graph of f(x) over [17, x1], [x1, x2], and [x2, 27.7] are all equal to 1/3 of the length
L found in Problem 2b.
Present numerical results in the following table:
x1
x2
Also, plot the points (xi, f(xi)), i = 1, 2 marked by ∗ together with the curve of
for the top profile of the noble beast. Again, make sure to scale the plot by
axis([0,30,0,30])
3. Halley’s comet last reached perihelion (closest to the sun) on Fe
uary 9th, 1986. Its
position and velocity components at this time were xy
z
=
XXXXXXXXXX− XXXXXXXXXX
XXXXXXXXXX
,
x′y′
z′
=
− XXXXXXXXXX− XXXXXXXXXX
− XXXXXXXXXX
where the position is measured in astronomical units (the earth’s mean distance from
the sun), and the time in years. The equations of motion for the comet are
x′′(t) = − µx(t)
(r(t))3
,
y′′(t) = − µy(t)
(r(t))3
,
z′′(t) = − µz(t)
(r(t))3
,
3
where
(t) =
√
(x(t))2 + (y(t))2 + (z(t))2, µ = 4π2,
and r(t) is distance between the comet and the sun.
(a) Solve the IVP of the first order system found above in the interval t ∈ [0, 100] (t
= 0 was the year of 1986) by the adaptive method (ode45 from Matlab). Set up
the parameters for the ode45 ODE solver as follows
options = odeset(’RelTol’,1e-8,’AbsTol’,1e-8);
Present numerical results in the following table:
x(100)
y′(100)
z′′(100)
(b) i. Use the numerical solution produced in Problem 3a and the cubic spline inter-
polation with the clamped boundary condition to find the position of Halley’s
comet at t∗ when the next closest approach of the comet to the sun. Present
your script for computing t∗ and present the numerical results in the following
table:
x(t∗)
y(t∗)
z(t∗)
ii. Use Matlab’s plot3 program to plot the trajectory of the comet together with
the position of the sun (x = 0, y = 0, z = 0) and the position of the comet at
the time t∗ found above.
Instructions for this plot:
• Mark the sun by a red ∗ and mark the position of the comet at the time
t∗ by a blue ∗.
• Adjust the plot by
view(-97.92, 23.64)
(c) Then compute the area between the curve of y(t) and the curve of z(t) over the
interval t ∈ [70, 80]. Use the numerical solution generated in Problem 3a and
the cubic spline interpolation with the clamped boundary condition to represent
functions y(t) and z(t) in this task.
4. This problem is to solve the following BVP (Boundary Value Problem) of a 4th orde
differential equation by the shooting method:
u′′′′(x) = f(x, u(x), u′(x), u′′(x), u′′′(x)), x ∈ (a, b), (1)
u(a) = b1, u
′(a) = b2, u(b) = b3, u
′(b) = b4. (2)
4
(a) Consider the related IVP (Initial Value Problem):
u′′′′(x) = f(x, u(x), u′(x), u′′(x), u′′′(x)), x ∈ (a, b], (3)
u(a) = b1, u
′(a) = b2, u
′′(a) = y1, u
′′′(a) = y2. (4)
For the following specific data:
f(x, u(x), u′(x), u′′(x), u′′′(x))
= 32cos(x)sin(x)− 2cos(cos(2x))sin(cos(2x)) + sin(u′(x)),
a = 1, b = 2, b1 = 0.9, b2 = −0.8, y1 = −3.6, y2 = 3.3,
use ode_rk33.m with Nh = 100 to solve the IVP described by equations (3) and
(4). Present the numerical solution in a table as follows:
x XXXXXXXXXX
u(x)
u′(x)
u′′(x)
u′′′(x)
(b) Consider the following function:
F(y) =
[
F1(y1, y2)
F2(y1, y2)
]
=
[
u(b)− b3
u′(b)− b4
]
, with y =
[
y1
y2
]
,
where u(x) is the solution to the IVP described by equations (3) and (4), b3, b4 are
specified in the boundary condition of the BVP described by equations (1) and
(2). Implement this function in Matlab with the following interface:
function F = two_point_ShMeth_beam_nonlinearF(y, f_fun, a, b, ...
1, b2, b3, b4, Nh, varargin)
In this Matlab function, use ode_rk33.m to solve for u(x) from the IVP described
y equations (3) and (4). The inputs
f_fun, a, b, b1, b2, b3, b4
are from the BVP, Nh is for the ODE solver ode_rk33.m.
Then, use the BVP described by equations (1) and (2) to test the Matla
two_point_ShMeth_beam_nonlinearF.m
with
f(x, u(x), u′(x), u′′(x), u′′′(x))
= 32cos(x)sin(x)− 2cos(cos(2x))sin(cos(2x)) + sin(u′(x)),
a = 1, b = 2, b1 = 0.1, b2 = 0.15, b3 = 0.2, b4 = 0.25.
5
Fo
y =
[
0.1
0.3
]
and Nh = 100, this Matlab function can produce
F(y) =
[
6.281888141405834e− 01
1.849268826683720e+ 00
]
The acceptable answer to this problem, i.e. Problem 4b should contain the Matla
function two_point_ShMeth_beam_nonlinearF.m, the script to compute F(y),
and a statement such as
This Matlab function of mine reproduces F(y) given above
(c) Use the shotting method to solve the BVP described by equations (1) and (2) with
f(x, u(x), u′(x), u′′(x), u′′′(x)) = 32cos(x)sin(x)
−2cos(u′′(x))sin(u(x)) + sin(u′(x)),
a = 1, b = 2, b1 = 0.1, b2 = 0.15, b3 = 0.2, b4 = 0.25.
Note that function f(x, u(x), u′(x), u′′(x), u′′′(x)) is different from the one in the
previous subproblems.
Present your script for the related computations and present numerical solution
in the following table:
x XXXXXXXXXX
u(x)
u′(x)
u′′(x)
u′′′(x)
5. Consider the following BVP for the steady state Burgers equation:
−y′′(x) + µy(x)y′(x) = f(x), x ∈ (0, 1),
y(0) = α, y(1) = β.
(a) Derive the finite difference equations for this BVP on a set of equispaced nodes:
a = x1 < x2 < · · · < xN = b, h = xi − xi−1, i = 2, 3, · · · , N.
Then describe the vector format F(u) = 0 for this system of finite difference
equations by identifying the component functions of F(u).
(b) Implement a Matlab function for F(u) and a Matlab function for JF(u). The
interfaces of these Maltab functions should be as follows:
6
function F = two_point_nl_F(u, x, f_fun, alpha, beta)
function J = two_point_nl_FJ(u, x, f_fun, alpha, beta)
(c) Write a Matlab script that uses both the Newton method and Broyden method
to solve the finite difference equation for the Burgers equation with the following
data:
f(x) = −ex−2x2
(
xex
(
− 1− x+ 2x2
)
+ ex
2(
2− 5x− 4x2 + 4x3
))
,
µ = 1, α = 0, β = 1, h = 0.002.
Present numerical results in the following table:
Method Number of iterations y(0.45)
Broyden
Newton
Some instructions for these computations:
• Use tol = 10^(-8) and nmax = 100 for stopping the iterations.
• Use the approximate Jacobian by the complex variable method for the initial
matrix B0 required by the Broyden’s method.
• Try the linear approximation as the initial guess for solving F(u) = 0.
(d) Write a Matlab script that uses the numerical solution generated in Problem 5c
and a suitable interpolation method to find an approximation to y(2/π). Jus-
tify your choice of the interpolation method and the related data for computing
this approximation to y(2/π) from the point of view of both the accuracy and
efficiency.
7