代写IVP留学生、Matlab编程语言调试、Matlab代写代做、system代写

日期：2019-12-15 10:29

4. Halley’s comet last reached perihelion (closest to the sun) on February 9th, 1986. Its

position and velocity components at this time were

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,

and r(t) is distance between the comet and the sun.

(a) (10 points) Reduce this IVP (initial value problem) of the system of second order

ODEs above to an IVP of a system of first order ODEs that can be solved by

numerical methods discussed in class.

(b) (10 points) 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:

(c) (15 points) Use the numerical solution produced in Problem 4b and cubic spline

interpolation with the clamped boundary condition to predict the position of Halley’s

comet at t = 55. Present the numerical results in the following table:

Also, 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 = 55 found above.

Instructions for this plot:

• Mark the sun by a red ∗ and mark the position of the comet at the time t = 55

by a blue ∗.

• Adjust the plot by

view(-116.9, 20.8)

(d) (10 points) Compute the length for comet’s trajectory plotted in Problem 4c.

(e) (20 points)

i. Use the numerical solution produced in Problem 4b to make a plot of the

function r(t), 0 ≤ t ≤ 100.

ii. Then use this plot and the numerical solution produced in Problem 4b to

esptimate/predict the time t

∗

, i.e., which year when the next closest approach

of the comet to the sun.

iii. Use the numerical solution produced in Problem 4b and the cubic spline interpolation

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:

iv. 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)

(f) (15 points) Locally use the numerical solution generated in Problem 4b and the

cubic Hermite interpolation to predict the position of Halley’s comet at t = 50.

Present your script for related computations and present the numerical results in

the following table:

It is required to describe and present the data used to construct the Hermite

interpolation.

6

(g) (10 points) Make a plot of y(t) and z(t) over the interval [70, 80]. Use the numerical

solution generated in Problem 4b and the cubic spline interpolation with

the clamped boundary condition to represent functions y(t) and z(t) in this task.

Scale the plot by

axis([70, 80, -4, 16])

(h) (10 points) Find the times in the time interval [70, 80] where the curves for y(t) and

z(t) intersect with each other. Use the numerical solution generated in Problem

4b and the cubic spline interpolation with the clamped boundary condition to

represent functions y(t) and z(t) in this task. Present your script for related

computations and present your solution in the following table:

t1

Add more rows if more than one such a time t are found.

(i) (10 points) Then compute the area between the curve of y(t) and the curve of z(t)

over the interval [70, 80].

5. This problem is to solve the following BVP (Boundary Value Problem) of a 4th order

differential equation by the shooting method:

u

0000(x) = f(x, u(x), u0

(x), u00(x), u000(x)), x ∈ (a, b), (1)

u(a) = b1, u0

(a) = b2, u(b) = b3, u0

(b) = b4. (2)

(a) (10 points) Consider the related IVP (Initial Value Problem):

u

0000(x) = f(x, u(x), u0

(x), u00(x), u000(x)), x ∈ (a, b], (3)

u(a) = b1, u0

(a) = b2, u00(a) = y1, u000(a) = y2. (4)

Reduce this IVP of a 4th order differential equation to the following IVP of first

order differential equations:

w0(x) = G(x, w), x ∈ (a, b], (5)

w(a) = w0. (6)

The formulas for G(x, w) and w0 should be presented as part of the solution.

(b) (10 points) Use ode_rk33.m with Nh = 100 to solve the IVP described by equations

(3) and (4) with the following information:

f(x, u(x), u0(x), u00(x), u000(x))= 32cos(x)sin(x) − 2cos(cos(2x))sin(cos(2x)) + sin(u0(x)),

a = 1, b = 2, b1 = 0.9, b2 = −0.8, y1 = −3.6, y2 = 3.3.

Present numerical solution in a table as follows:

(c) (10 points) Consider the following function:

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, ...

b1, b2, b3, b4, Nh, varargin)

In this Matlab function, use ode_rk33.m to solve for u(x) from the IVP discussed

in Problem 5a. 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 Matlab

two_point_ShMeth_beam_nonlinearF.m

with

f(x, u(x), u0(x), u00(x), u000(x))= 32cos(x)sin(x) − 2cos(cos(2x))sin(cos(2x)) + sin(u0(x)),

a = 1, b = 2, b1 = 0.1, b2 = 0.15, b3 = 0.2, b4 = 0.25.

and Nh = 100, this Matlab function can produce

6.281888141405834e − 01

1.849268826683720e + 00

The acceptable answer to this problem, i.e. Problem 5c should contain the Matlab

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

8

(d) (20 points) Solve the BVP described by equations (1) and (2) with

f(x, u(x), u0(x), u00(x), u000(x)) = 32cos(x)sin(x)−2cos(u00(x))sin(u(x)) + sin(u0(x)),

a = 1, b = 2, b1 = 0.1, b2 = 0.15, b3 = 0.2, b4 = 0.25.

Note that function f(x, u(x), u0

(x), u00(x), u000(x)) is different from the one in the

previous subproblems.

This can be done by 2 steps:

Step 1: Compute the zero of F(y) by the Broyden method. In this computation, use

approxJ_compl.m to compute the initial matrix B0 with

y0 = [0.1; 0.3]; EPS = 10^(-6)

and use Nh = 100 for the ODE solver used in

two_point_ShMeth_beam_nonlinerF.m

Then use the following inputs for the Broyden method:

y0 = [0.1; 0.3]; tol = 10^(-12); nmax = 200;

Step 2: Then, use ode_rk33.m with Nh = 100 to solve the IVP described by (3) and

(4) in which y1 and y2 are found in Step 1.

Present your script for the related computations and present numerical solution

in the following table:

6. Consider the following BVP for the steady state Burgers equation:

(a) (5 points) 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) (5 points) Derive the Jacobian JF(u) for the nonlinear function F(u).

(c) (10 points) Implement a Matlab function for F(u) and a Matlab function for JF(u).

The interfaces of these Maltab functions should be as follows:

9

function F = two_point_nl_F(u, x, f_fun, alpha, beta)

function J = two_point_nl_FJ(u, x, f_fun, alpha, beta)

(d) (10 points) 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:

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.

(e) (5 points) Write a Matlab script that uses the numerical solution generated in Problem

6d and a suitable interpolation method to find an approximation to y(1/π).

Justify your choice of the interpolation method and the related data for computing

this approximation to y(1/π) from the point of view of both the accuracy and

efficiency.

7. Consider the following IBVP for the Burgers equation:

∂y(x, t)∂t −∂2y(x, t)

∂x2 + q(x, y(x, t))∂y(x, t)

∂x = f(x, t, y(x, t)),

x ∈ (0, 1), t ∈ (0, 5],

y(x, 0) = y0(x), x ∈ (0, 1),

y(a, t) = α(t), y(b, t) = β(t), t ∈ (0, 1].

with

q(x, t, y) = y,

f(x, t, y) = 3cos(t + x(1 − x)) − (2x − 1)(1 − 2x + cos(t + x(1 − x)))y,

y0(x) = sin(x − x2), α(t) = sin(t), β(t) = sin(t),

(a) (5 points) Derive the ODE system based on using finite difference to approximate

the x partial derivatives for this IBVP on a set of equispaced space variable nodes:

a = x0 < x1 < x2 < · · · < xNx < xNx+1 = b,

h = xi − xi−1, i = 1, 3, · · · , Nx + 1.

Then put these ODEs together in the form: u

0

(t) = F(t, u).

10

(b) (5 points) Follow the Matlab function heat_eq_nonlinear_ode_sys in lecture

slides to implement a Matlab function for the ODE system in Problem 7a such

that it has the following interface:

function F = heat_eq_nonlinear_ode_sys_Burgers(t, u, ...

x, q_fun, f_fun, alpha_fun, beta_fun)

(c) (10 points) Write a Matlab script that uses the Crank-Nilcon ODE solver ode_CN_Broyden

together with the function heat_eq_nonlinear_ode_sys_Burgers implemented

above to solve this IBVP with h = 0.01. The input variable Nt for the CrankNicolson

ODE solver should be chosen such that the step size τ for the time variable

t is the same as h = 0.01. Present your numerical solution in the following table

y(0.5, 2.25) ∂y(0.5, 2.25)

∂t

(d) (5 points) Write a Matlab script that uses the numerical solution generated above

to find an approximation to y(2/π, π).

(e) (5 points) Write a Matlab script that makes a surface plot of y(x, t) for (x, t) ∈[0, 1] × [0, 5]. The script should label both the x axis and the t axis, and it should

use the following Matlab command to adjust the viewing angle of this surface plot:

view(-1.1e02, 1.9e01)

(f) (5 points) Write a Matlab script for making a movie that displays the curves for

y(x, t) for a chosen sequence of values of t.

11