The initial value problem for ordinary differential equations of the previous labs is only one of the two major types of problem for ordinary differential equations. The other type is known as the ``boundary value problem'' (BVP). A simple example of such a problem would describe the shape of a rope hanging between two posts. We know the position of the endpoints, and we have a second order differential equation describing the shape. If the two conditions were both given at the left endpoint, we'd know what to do right away. But how do we handle this ``slight'' variation?
This lab is concerned with two of the most common approaches to solving BVPs as well as a combined IVP-BVP for a partial differential equation. The extra credit problem introduces a third approach to solving BVPs. The discussion in this lab is limited to relatively simple approaches in a single space dimension and is intended to give the flavor of these approaches, each of which could easily be the subject of a full semester's course. Except for the extra credit exercise, these methods are easily extended to two and three space dimensions.
The approaches included in this lab are the following:
A one-dimensional boundary value problem (BVP), is similar to an initial value problem, except that the data we are given isn't conveniently located at a starting point, but rather some is specified at the left end point and some at the right. (We're also usually thinking of the independent variable as representing space, rather than time, in this setting).
We will be using the following problem as the illustrative example
for several of the following exercises.
The clothesline BVP: A rope is stretched between two
points. If the rope were weightless, or if it were rigid, it would lie along a
straight line; however, the rope has a weight and is elastic, so it
sags down slightly from its ideal linear shape. We wish to determine the curve
described by the rope. We will use the variable 
to denote
horizontal distance and 
to denote height of the rope
at the point
.
The equation for the curve described by the rope can be derived (this
is not a proof: it is a description of why you should believe the
equation) in the following manner. Suppose that the tension in the
rope is 
(a constant because the rope is in equilibrium), and
consider a tiny piece of the rope of length 
, and
with mass per unit length 
. The total mass of the differential
piece of rope is 
, so that the force due to gravity is
directed downward and is given by
. This piece of rope observes
forces on each of its ends. The magnitudes of these forces are equal
to the tension,
, and the directions are given by the slope of the
curve at the ends of the differential piece. Hooke's law says that
the tension is proportional to the amount of strain in the string,
, where 
is a constant of proportionality (Young's modulus).
Hence, the equation can
be written as
![$\displaystyle -K\left.\frac{du}{dx}\right]_{\mbox{right}}+K\left.\frac{du}{dx}\right]_{\mbox{left}}
=- \rho g dx. $](img11.png)
Dividing both sides by
and letting
yields the equation
. There is no reason that the ``constant'' of proportionality cannot change from place to place. For the sake of definiteness, assume
varies as 
, for constant 
, representing a rope (or spring)
whose stiffness varies from end to end. The result is the equation
When
, this equation is called the ``Poisson equation'' and also describes
the distribution of heat in a solid bar, among other common physical problems.
For the sake of definiteness, take
and 
, the left end of
height 
at 
, and the right end height of 
at 
.
Thus, the system to be solved is
The ``derivation'' presented above for the shape of the
rope is suggestive of a way to solve for the shape, called
the ``finite difference method.'' Assume that
we have divided the interval up into 
equal intervals of
width 
determined by 
points. Denote the spatial
points 
,
. Approximate the
value of 
by 
. Also approximate the Young's modulus
function as
.
Now, consider the
interval as if it were the differential
piece of rope mentioned in the derivation. Using the standard finite difference
approximation for a derivative, the slope of the rope at the left of
the
interval could be approximated as
and the slope on the right of the
interval could be
approximated as
. The difference between
these is an approximation of the second derivative
Along similar lines, approximate the first derivative as
Both approximations (of
and 
) have the same
Taylor-series (truncation) error of
.
Put all these together into (1) to get
and, as above,
and
.
We can associate this equation with the solution value at
,
except for 
and 
(do you see why?). Conveniently, those
are the points at which we have boundary conditions specified.
In particular, let us look at approximating our rope BVP at 6
points. We set up the ODE at points 1, 2, 3, and 4, and associate the
boundary conditions with the
and 
solution values. Note
that
. I also
multiplied through by
to make things look nicer:
Actually, in Equation (2), the quantities 
and 
are not really variables, being fixed by the boundary conditions. Hence
the only variables are 
, 
, 
and 
. The system can
be rewritten as
and this system has been formatted to suggest the matrix equation
By discretizing the differential equations we have created a set of
linear algebraic equations that have the symbolic form 
.
To set up and solve the equations
(3) in Matlab, we could type:
![$\displaystyle \left[\begin{array}{rrrr} -2(1+c\Delta x)&+(1+1.5c\Delta x)&+0&+0...
...Delta x^2 0.4\Delta x^2 0.4\Delta x^2-(1+4.5c\Delta x)u_5\end{array}\right]$](img70.png)
N = 4;
C = 0.05;
RHOG = 0.4;
% N interior mesh points, N+1 intervals
dx = 5.0 / ( N + 1 );
x = dx * (0:N+1);
A = [ -2*(1+C*dx) +(1+1.5*C*dx) 0 0;
+(1+1.5*C*dx) -2*(1+2*C*dx) +(1+2.5*C*dx) 0;
0 +(1+2.5*C*dx) -2*(1+3*C*dx) (1+3.5*C*dx);
0 0 +(1+3.5*C*dx) -2*(1+4*C*dx) ];
ULeft=1;
URight=1.5;
b = [ RHOG*dx^2-(1+0.5*C*dx)*ULeft
RHOG*dx^2
RHOG*dx^2
RHOG*dx^2-(1+4.5*C*dx)*URight];
U = A \ b;
U = [ULeft; U; URight]
Make sure you understand the first and last components in b.
You should recall that the backslash notation
is shorthand for saying U=inv(A)*b
but tells Matlab to solve the equation A*U=b
without actually forming the inverse of A.
Remark: The vector x is a row vector and the vector U is a column vector! This is the convention that has been followed for the _ode.m files and will be followed throughout these labs.
,
,
into your chosen
equation. Show the result is essentially zero.
is upper-case C,
is dx and
is U(1:6)
in exer1a.
would be a solution if 
and 
.
As a second verification step, make a copy of exer1a.m
called exer1c.m with rhog=0 and URight=1
and check that the discrete solution is (exactly or to roundoff) correct.
would be a solution if 
, 
and 
.
As a third verification step, make a copy of exer1a.m
called exer1d.m with rhog=C and and with
boundary values ULeft and URight chosen to match
the solution you are testing.
Check that the discrete solution is (exactly or to roundoff) correct.
would be a solution if
, 
and 
.
As a fourth verification step, make a copy of exer1a.m
called exer1e.m with constant rhog, which
is used four times, replaced by the four values of the
vector 2+4*c*x(2:5)' and with
boundary values chosen to agree with
and check that the discrete solution is (exactly or to roundoff) correct.
and 
on its ends. Please include this plot with your summary.
Remark: The exact solutions 
, 
, and 
can be used as exact discrete solutions and verification tests
only when the approximation expressions for first and second
derivatives are sufficiently accurate. Because the mesh is uniform,
the approximate expressions used here for the derivatives yield
the same values as using the usual continuous expressions so long as
is a quadratic (or lower) polynomial, so there is no truncation
error and solutions are correct to roundoff.
The purpose of the previous exercise is to verify that you have copied the code correctly and to illustrate the powerful verification strategy of checking against known exact discrete solutions. You should always use a small, simple problem to verify code by comparison with hand calculations. If possible, you should also compare results with theoretical results and with results achieved using a different method. In the next exercise you will be modifying the above code to handle the case of large N and solving a slightly more realistic problem. Of course, you don't want to bother typing in the matrix A if N is 100 or more, so you will be writing Matlab code to do it.
function [x,U] = rope_bvp(N) % [x,U] = rope_bvp(N) % comments % your name and the date
as you did for
exer1c.m above. Check your results for N=4.
Fix any mistakes before continuing.
and check the terms against your code, one term at a time.
as you did for
exer1d.m above. Check your results for N=4.
Fix any mistakes before continuing.
and check the terms against your code, one term at a time.
as you did for
exer1d.m above. Check your results for N=4.
Fix any mistakes before continuing.
and check the terms against your code, one term at a time.
Remark: You may wonder why the four-point mesh solution
and the 119-point mesh solution seem to agree at the four common points.
This behavior is highly unusual. It happens because the
solution of the differential equation is almost quadratic and the
difference scheme exactly reproduces quadratic functions. For larger
values of
, the solution looks less like a quadratic, and the
solution for 
agrees less well with the solution for 
.
Try it, if you wish.
In the previous section, you saw an example of the finite difference
method of discretizing a boundary value problem. This method is based
on a finite difference expression for the derivatives that appear in
the equation itself. The finite difference method results in a
list of values that approximate the true solution at the set of mesh
points. Approximate values between the mesh points might be generated
using interpolation ideas, but the method itself does not depend on
any such interpolation. The reason that
was chosen for comparison
with
in the previous exercise is because the four
values in the
case appear among the 
-values in the
case.
An alternative approach, called the ``finite element method'' (FEM) is based on approximating the unknown as a sum of simple ``shape functions'' defined over the mesh intervals. Since the finite element solution is actually a function, it is defined over the same spatial interval as the true solution and much of the machinery of functional analysis is available for proving facts about the method and solutions that arise. As a consequence, the FEM occupies a large part of the mathematics literature. You can find the FEM discussed in Quarteroni, Sacco, and Saleri, Sections 12.4 and 12.5.
In this section, you will see the FEM applied to a particular boundary value problem. The problem is somewhat simpler than the clothesline problem discussed above, but contains the same essential features. Consider the equation
in the interval ![$ [0,1]$](img97.png){kind=small}
, for 
a given function,
and with boundary values
, assumed to satisfy the same boundary conditions
(5), and integrating some terms by parts. In this case,
the weak form is given by
![$\displaystyle -\int_0^1 y'(x)v'(x)dx + \left[\rule{0pt}{12pt}y'(x)v(x)\right]_0^1 +\int_0^1 y'(x)v(x)dx +\int_0^1y(x)v(x)dx=\int_0^1f(x)v(x)dx$](img101.png) |
satisfies (5), the bracketed term
drops out and the result is
To approximate the function 
, choose an odd integer
and a
set of functions 
for
, defined on the interval
, that form a basis
of some reasonable approximating function space. For most finite element
constructions, these functions satisfy the following characteristics:
subintervals, each of length 
,
so that a sequence of spatial points is given by 
for
(
and 
). The quadratic
Lagrange functions are defined in the following way
(see Quarteroni, Sacco, Saleri, p. 562).
and it has good approximation properties. It is also true
that each of the functions
, for
satisfies
the boundary conditions (5).
Assume that an approximate solution to (6) can be written as
.
Plugging (9) into (6)
and choosing
yields
equations of the form
as the components of a (column) vector
, the values 
as the components of a matrix
, and the values 
as the components of a (column)
vector
, then (10) can be written as the
matrix equation
In the following exercises, you will write Matlab functions to construct
the basis functions
in (7) and (8),
to evaluate the matrix elements
and vector components
in (10), and solve the matrix equation
(11).
Remark: When the FEM is programmed, the integrals in (11) are typically performed element-by-element. This is particularly important in multidimensional cases. Nontheless, we will be using a conceptually simpler approach to the integrations.
used in (7) and (8) will be
evaluated as 
, where
. The expression 
is valid
even when 
, although Matlab does not allow subscripts equal to zero.
In the Matlab functions below, you should regard the variable x as a scalar value, not a vector. Attempting to write vector (componentwise) code only complicates matters here.
by completing the
following outline
![\begin{displaymath}\left\{
\begin{array}{ll}
\frac{(x-x_{n-1})(x-x_{n-2})}{(x_n-...
...x \leq x_{n+2} [5pt]
0 & \text{ otherwise}
\end{array}\right.\end{displaymath}](img113.png)
![$\displaystyle \left\{\begin{array}{ll}
\frac{(x_{n+1}-x)(x-x_{n-1})}{(x_{n+1}-x...
...or } x_{n-1}\leq x\leq x_{n+1} [5pt]
0 & \text{ otherwise}
\end{array}\right.$](img114.png)
function z=phi(n,h,x)
% z=phi(n,h,x)
% Lagrange quadratic basis functions
% your name and the date
if numel(x) > 1
error('x is a scalar, not a vector, in phi.m');
end
if mod(n,2)==0 % n is even
if (n-2)*h < x & x <= n*h
z= ??? code implementing first part of (7) ???
elseif n*h < x & x <= (n+2)*h
z= ??? code implementing second part of (7) ???
else
z=0;
end
else % n is odd
??? code implementing (8) ???
end
N=7; h=1/(N+1); x=linspace(0,1,97); mesh=linspace(0,1,N+2); for k=1:numel(x) y3(k)=phi(3,h,x(k)); y4(k)=phi(4,h,x(k)); end plot(x,y3,'b') hold on plot(x,y4,'r') plot(mesh,zeros(size(mesh)),'*') hold offYou should observe that each
takes the value 1 at a single
mesh node (x(k), indicated with an asterisk), takes the value
zero at all other mesh nodes, is continuous, and is parabolic or zero
between any two mesh nodes. Please include this plot
with your summary file.
is a continuous function by showing that
the pieces match up at 
, 
, and 
for even 
and at 
and 
for odd
. (You don't need Matlab
to do this.) Similarly, show that
 |
. Differentiate (7) and (8)
by hand to find
and use your
formulæ for the function phip.m with signature
function z=phip(n,h,x) % z=phip(n,h,x) % derivative of Lagrange quadratic basis functions % your name and the date
and
(multiply
by 
to get a better scaling) on the
same plot. Examine the plot carefully and convince yourself that
appears to be the derivative of 
.
and 
. Examine both cases from
(7).
, use the finite difference
expression
 |
to estimate
. Does it agree up to
roundoff with the result from phip? Similarly for
.
In the following three exercises, you will write m-files to construct and verify the three pieces of the matrix A in (10). Following that, you will construct the full matrix A and solve for the finite element solution of the given problem (4).
,
In order to do the integrations, you will be using code that I give you that provides a uniform way to do the required integrations.
function q=gaussquad(f1,k1,f2,k2,h,a,b)This gaussquad function will provide the exact value, not merely an approximate value, of the integral for the cases considered in this lab: piecewise low degree polynomials.
in (12). Call the resulting
matrix A1. Do not overlook the fact that the basis function
derivatives appear in (12), not the basis functions themselves, and
there is that pesky minus sign in front of the integral. Please
include the values of A1 in your summary file.
is contained in the
interval
![$ [(n-2)h,(n+2)h]$](img153.png){kind=small}
. You could use this fact to shrink
your integration limits and improve the efficiency of the integration,
but it is not required.
and of 
do not intersect when
, you can take advantage of this fact to
avoid computing components
that must be zero, but it is
not required.
norm(A1-A1','fro')is zero or roundoff. Add this code to exer4.m.
for the 97
values x=linspace(0,1,97) and plot it. You should observe that
it is equal to 1 except near the endpoints of the interval. As a consequence,
it has zero derivative, except near the endpoints of the interval. Please
include the plot with your summary.
, and you just saw that 
is zero
except for
and
,
verify the zeros in positions 
.
, you
should be able to see why the following code
N=7; h=1/(N+1); v=(1:N)'*h; % v=x A1*vshould yield a vector that is zero except in the positions
.
Add this code to exer4.m and verify that it does.
with
has solution
. You can solve this BVP using FEM.
Write a Matlab function m-file with the signature
function z=rhs4(n,h,x) % z=rhs4(n,h,x) % your name and the dateto compute the constant function equal to (-2) everywhere. (This is an almost trivial exercise. It is needed so that gaussquad.m can be used.) Add code to exer4.m to use the gaussquad.m function to compute the vector components
and call the resulting vector
RHS4. Since the quadratic function
satisfies the
boundary conditions, it can be written exactly as a combination
of the 
! Hence, the following code should yield the zero vector.
N=7; h=1/(N+1); xx=(1:N)'*h; % the variable x has already been used. v=xx.*(1-xx); v-A1\RHS4 % should be zeroAdd this code to exer4.m. Please include the values of RHS4 in your summary.
The tests in exer4.m construct and test the matrix A1, so you should now be reasonably sure that A1 is correct. In the following exercise, you will construct and test A3. In the subsequent exercise, you will construct and test A2 so that you have all of the matrix A.
 |
has the exact solution 
, and this quadratic
function can be expressed exactly as a sum of the
.
function z=rhs5(k,h,x) % z=rhs5(k,h,x) % your name and the dateto compute the function
. Add code to exer5.m
to use rhs5.m and gaussquad.m
to compute the (column) vector
, and call
the resulting vector RHS5. Include these values in your summary
at each of the nodes
exactly, with only roundoff errors.
Add code to exer5.m to check that this is true.
If it is not true, there is a mistake somewhere. Fix it before continuing.
has exact solution
.
function z=rhs6(k,h,x) % z=rhs6(k,h,x) % your name and the dateto compute the function
.
Add code to exer6.m
to compute the (column) vector
, and call
the resulting vector RHS6.
exactly, with only roundoff errors.
Remark: Writing test scripts like exer4.m, exer5.m and exer6.m allows you to re-test your work at any time. This strategem helps you maintain confidence in your code and also allows you to propose and test modifications easily.
 |
has exact solution
, with
.
.
function [x,Y]=solve7(N) % [x,Y]=solve7(N) % ... more comments ... % your name and the datethat performs the following tasks:
, and
take the ratio as the error(h) divided by error(h/2).
N h error ratio 7 1.2500e-1 ________ ________ 15 6.2500e-2 ________ ________ 31 3.1250e-2 ________ ________ 61 1.6129e-2 ________ ________ 121 8.1967e-3 ________Remark: The rate of convergence appears higher than expected from theory. This rate is an artifact of the way the error is computed. To do it properly, one needs to compute
integral
errors that involve more than just the node points. Convergence at
the node points is one higher order than expected, a feature called
``superconvergence.''
A partial differential equation (PDE) involves
derivatives of a function
which depends on more than one
independent variable. One interesting PDE is the one-dimensional
Burgers' equation, which is a one-dimensional nonlinear equation
whose nonlinear term is similar to the one in the Navier-Stokes equations
of fluid flow. In this case, the variable will be called
and it
is a function of space (
) and time (
), 
. The variable
is called ``velocity'' in the Navier-Stokes equations.
Burgers' equation on the spatial interval
can be written as
is a constant. Boundary conditions can be taken as
and 
, both for all time. Two space boundary conditions
are necessary because the equation is second-order in space.
Since the equation is first-order in time, only one intitial condition
is needed.
To get an idea of what the solution to Burgers' equation might look like,
first imagine that 
and also that the coefficient 
is a constant,
so that the equation becomes the wave equation
If
, this equation represents a right-going wave that moves without
changing shape. To see this, note that for any given function of one
variable, 
, 
is a solution of the wave equation. Changing
to a small, positive number means that the wave propagates as before,
but slowly spreads and decays to zero. Finally, since the coefficient 
is
not constant, the wave propagates faster where
is larger and
slower where
is smaller. Thus a wave that is larger to the left and
smaller to the right will steepen as it propagates to the right.
Look at (16) and pretend, just for a moment, that
time is frozen. The equation suddenly starts looking like the
boundary value problem (1) that we solved before,
only with a extra
term on the left and with
replacing the right
side. This observation is the
basis for the ``method of lines,'' wherein the spatial
discretization is performed separately from the temporal. Consider
the function 
, but think of it as a
function of
first. The resulting BVP can be
written with primes denoting spatial differentiation so it
looks more like what we have been doing.
where we are focussing on the left side for the moment, along with spatial boundary conditions
and 
We solved a BVP a lot like this one above. We broke the interval
into 
subintervals and labelled the
resulting points
. Next, we
defined 
as being the approximate
solution at
. Keep in the back of you mind,
though, that
is really still a function of
,
. We will denote the vector of values
by 
, because we
have already used the unsubscripted
to denote the continuous
solution. Keep in mind that 
is a function of time. We will
use the same finite difference discretization as before for the term 
.) The
resulting discrete equations become
But we remain on familiar ground: (19) is just a system of IVPs!
We know a bunch of different methods to solve it. It turns out that
the system is moderately stiff, and will become more stiff when N
is taken larger and larger. We will use backwards Euler to solve this
system. Recall that backwards Euler requires an m-file to evaluate
both the function
and its partial derivative (Jacobian)
In Matlab notation, the variable
will be a matrix whose
entries 
approximate the values 
. In Matlab
notation, for the
time interval, U(:,k)
represents the (column vector of) values at the
locations x(:). The initial condition for U
should be a column vector whose values are specified at
the locations x(:).
and using N=500.
refers to the right side of (19).
Here is an outline:
function [F,D]=burgers_ode(t,U)
% [F,D]=burgers_ode(t,U)
% compute the right side of the time-dependent ODE arising from
% a method of lines reduction of Burgers' equation
% and its derivative, D.
% Boundary conditions are fixed =0 at the
% endpoints x=0 and x=1.
% A fixed number of spatial points (=N) is used.
% The variable t is not used, but is kept as a place holder.
% U is the vector of the approximate solution at all spatial points
% output F is the time derivative of U (column vector)
% output D is the Jacobian matrix of
% partial derivatives of F with respect to U
% spatial intervals
N=500;
NU=0.001;
dx=1/(N+1);
ULeft=1; % left boundary value
URight=0; % right boundary value
F=zeros(N,1); % force F to be a column vector
D=zeros(N,N); % matrix of partial derivatives
% construct F and D in a loop
for n=1:N
if n==1 % left boundary
F(n) = ??? Function, left endpoint ???
D(?,?)= ??? Derivative, left endpoint ???
elseif n<N % interior of interval
F(n) = ??? Function, interior points ???
D(?,?)= ??? Derivative, interior points ???
else % right boundary
F(n) = ??? Function, right endpoint???
D(?,?)= ??? Derivative, right endpoint???
end
end
F(n)= ??? Function, interior points ???
with the discretization for F(n) given in (19).
D(?,?)= ??? Derivative, interior points ???
with the values D(n,n), D(n,n+1), and D(n,n-1)
according to the formula that was given above in (20)
and is repeated here.
 |
(21) |
will take on the three values n, n+1,
and n-1. To do, for example, D(n,n+1), write out
the expression for F(n) and differentiate it with respect
to U(n+1). Remark: It is easy to see from the formula that
for 
or 
.
F(n) = ??? Function, left endpoint ???
D(?,?)= ??? Derivative, left endpoint ???
with the expressions for F(n), D(n,n) and D(n,n+1).
F(n) = ??? Function, right endpoint???
D(?,?)= ??? Derivative, right endpoint???
with the expressions for F(n), D(n,n) and D(n,n-1).
N=500; % must agree with value inside burgers_ode.m x=linspace(0,1,N+2); x=x(2:N+1);and we will assume an initial velocity distribution that looks like a shallowly sloped wave.
UInit=((1-x).^3)';
results burgers_ode(0,UInit)
n F(n) D(n-1,n) D(n,n) D(n+1,n)
1 2.9761711475 -499.01396011 498.51296011
2 2.9465759259 1.99800798 -499.02590030 497.02789232
250 0.0976958603 219.12262501 -501.24899902 282.12637400
499 2.395210e-05 251.00094622 -502.00194821 251.00100199
500 1.197605e-05 251.00098406 -502.00198406
[t,U]=back_euler('burgers_ode',[0,1],UInit,100);
The solution should converge at each step. If you get the
``failed to converge'' message, your derivative (D)
is probably wrong. To help with grading, please include the value
of U(200,50) in your summary.
plot(x, U(:,k) )where x is the value assigned above. Please include this plot for the choice k=50 with your summary.
plot(x,U(:,1)) axis([0,1,0,1.5]) for k=2:100 pause(0.1); plot(x,U(:,k)); axis([0,1,0,1.5]) endYou should be able to see the ``wave'' steepen as it moves to the right. Please include the final frame of this sequence with your summary.
.
Shooting methods solve a BVP by reformulating it as an IVP, using initial values as parameters. The BVP is solved by finding the parameters that reproduces the desired boundary values. For this exercise, we return to considering the BVP of a hanging rope.
Taking as our example the rope BVP (1),
how much violence do we have to do to it
in order to make it look like an IVP? Well, we expect to have two
conditions at the left point and none at the right point. So let's
temporarily consider the related problem, where we have made up an
extra boundary condition at our initial value of
:
The following exercise attacks this system using a method called ``shooting.'' The strategy behind this method is the following.
, we can solve this problem, and get a
numerical solution at a sequence of points up to the right endpoint.
determines a (numerical) solution 
, we can regard the
difference between the value we got, and the value we want, as a
function
. (Where we write
,
but we should really write something like
, to emphasize that the
solution depends on the parameter.)
. The Matlab function
fzero will be used to find 
.
As you know, the ODE (1) can be written as a system
of first order differential equations. This is done by
identifying 
and 
, yielding the system
 |
 |
 |
|
 |
|
 |
For this exercise, you will be using built-in Matlab ODE routines to solve the initial value problem that arises from a guessed initial condition and then you will use the Matlab function fzero to solve for the correct initial condition.
function fValue=rope_ode(x,y) % fValue=rope_ode(x,y) computes the % rhs of the first-order system % your name and the dateNote that since we plan to use ode45 there is no need to add the Jacobian matrix to rope_ode.m! This is a great convenience, but in many cases you would have to provide a function for the Jacobian or ode45 (or ode15s, etc.) might fail. In that case, setting an option allows the Jacobian computation.
and use the Matlab
function ode45 to solve the system (22) on the
interval [0,5]. What is the value of 
? Is it positive
or negative? (Be careful: which of the components of y
corresponds with
?)
for which
the value of
is of the opposite sign as for
.
The correct value of
lies between the two values you just
found.
function F = rope_shoot ( alpha ) % F = rope_shoot ( alpha ) % comments % your name and the dateand this code should do the following:
;
![$ x\in[0,5]$](img242.png){kind=small}
;
.
, the function you just wrote will return
. When
is just right, it will return 0.
that makes
fzero
requires two parameters, a function handle (
@) first and
second the vector [alpha1,alpha2]
of the two values of
that you just found
for which 
has opposite signs. What is the value of
alpha you found?
and a height of 1.5 at
?
You do not need to send me this plot.
you just computed?
Back to MATH2071 page.