没有合适的资源?快使用搜索试试~ 我知道了~
lab04.pdf
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 112 浏览量
2023-06-18
13:13:35
上传
评论
收藏 118KB PDF 举报
温馨提示
试读
17页
lab04.pdf
资源推荐
资源详情
资源评论
MATH2071: LAB 4: BVPs and PDEs
Introduction Exercise 1
Boundary Value Problems Exercise 2
Discretizing a BVP Exercise 3
Finite element method Exercise 4
Burgers’ Equation Exercise 5
The Method of Lines Exercise 6
Extra Credit: Shooting Methods Exercise 7
Exercise 8
Extra Credit
1 Introduction
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 “b oundary value
problem” (BVP). A simple example of such a problem would describe the shape of a rope hanging between
two posts. We know the p osition 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:
• The Finite Difference method (FDM),
• The Finite Element method (FEM),
• The Method of Lines, and,
• The Shooting method (extra credit).
This lab will take four sessions. If you print this lab, you may prefer to use the pdf version.
2 Boundary Value Problems
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
1
from its ideal linear shape. We wish to determine the curve described by the rope. We will use the variable
x to denote horizontal distance and u(x) to denote height of the rope at the point x.
Forces on a clothesline
s
s
s
-
T
@
@
@I
T
?
ρg dx
dx
-
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
T (a constant because the rope is in equilibrium), and consider a tiny piece of the rope of length dx, and
with mass per unit length ρ. The total mass of the differential piece of rope is ρ dx, so that the force due to
gravity is directed downward and is given by ρg dx. This piece of rope observes forces on each of its ends.
The magnitudes of these forces are equal to the tension, T , 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 prop ortional to the amount
of strain in the string, T = −Kdu/dx, where K is a constant of proportionality (Young’s modulus). Hence,
the equation can be written as
−K
du
dx
right
+ K
du
dx
left
= −ρg dx.
Dividing both sides by dx and letting dx → 0 yields the equation
−
d
dx
K
du
dx
= −ρg
. There is no reason that the “constant” of proportionality cannot change from place to place. For the sake
of definiteness, assume K varies as K(x) = 1 + cx, for constant C, representing a rope (or spring) whose
stiffness varies from end to end. The result is the equation
(1 + cx)u
′′
+ cy
′
= ρg.
When c = 0, 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 ρg = 0.4 and c = 0.05, the left end of height 1 at x = 0, and the right
end height of 1.5 at x = 5. Thus, the system to be solved is
(1 + cx)u
′′
+ cu
′
= ρg
c = 0.05
ρg = 0.4 (1)
u(0) = 1
u(5) = 1.5
The ODE is linear. Linearity implies good things such as the existence and uniqueness of solutions.
2
3 Finite Difference Method for a BVP
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 N equal intervals of
width ∆x determined by N + 1 points. Denote the spatial points x
n
, n = 0, 1, . . . , N + 1. Approximate the
value of y(x
n
) by y
n
. Also approximate the Young’s modulus function as K
n
= K(x
n
) = (1 + cx
n
).
Now, consider the n
th
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
n
th
interval could be approximated as (y
n
−y
n−1
)/∆x and the slope on the right of the n
th
interval could be
approximated as (y
n+1
−y
n
)/∆x. The difference between these is an approximation of the second derivative
y
′′
n
≈
y
n+1
− 2y
n
+ y
n−1
∆x
2
.
Along similar lines, approximate the first derivative as
y
′
n
≈
y
n+1
− y
n−1
2∆x
.
Both approximations (of y
′
and y
′′
) have the same Taylor-series (truncation) error of O(∆x
2
).
Put all these together into (1) to get
(1 + cx
n
)
y
n+1
− 2y
n
+ y
n−1
∆x
2
+ c
y
n+1
− y
n−1
2∆x
= ρg
and, as above, c = 0.05 and ρg = 0.4. We can associate this equation with the solution value at y
n
, except
for n = 0 and n = N + 1 (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 n = 0 and n = 5 solution values. Note that
x
n
= n∆x. I also multiplied through by ∆x
2
to make things look nicer:
u
0
= 1
(1 + c∆x − c∆x/2)u
0
− 2(1 + c∆x)u
1
+ (1 + c∆x + c∆x/2)u
2
= 0.4∆x
2
(1 + 2c∆x − c∆x/2)u
1
− 2(1 + 2c∆x)u
2
+ (1 + 2c∆x + c∆x/2)u
3
= 0.4∆x
2
(2)
(1 + 3c∆x − c∆x/2)u
2
− 2(1 + 3c∆x)u
3
+ (1 + 3c∆x + c∆x/2)u
4
= 0.4∆x
2
(1 + 4c∆x − c∆x/2)u
3
− 2(1 + 4c∆x)u
4
+ (1 + 4c∆x + c∆x/2)u
5
= 0.4∆x
2
u
5
= 1.5
Actually, in Equation (2), the quantities u
0
and u
5
are not really variables, being fixed by the boundary
conditions. Hence the only variables are u
1
, u
2
, u
3
and u
4
. The system can be rewritten as
−2(1 + c∆x)u
1
+(1 + 1.5c∆x)u
2
+0 +0 = 0.4∆x
2
− (1 + 0.5c∆x)u
0
(1 + 1.5c∆x)u
1
−2(1 + 2c∆x)u
2
(1 + 2.5c∆x)u
3
+0 = 0.4∆x
2
0 +(1 + 2.5c∆x)u
2
−2(1 + 3c∆x)u
3
+(1 + 3.5c∆x)u
4
= 0.4∆x
2
0 +0 +(1 + 3.5c∆x)u
3
−2(1 + 4c∆x)u
4
= 0.4∆x
2
− (1 + 4.5c∆x)u
5
and this system has been formatted to suggest the matrix equation
−2(1 + c∆x) +(1 + 1.5c∆x) +0 +0
(1 + 1.5c∆x) −2(1 + 2c∆x) (1 + 2.5c∆x) +0
0 +(1 + 2.5c∆x) −2(1 + 3c∆x) +(1 + 3.5c∆x)
0 +0 +(1 + 3.5c∆x) −2(1 + 4c∆x)
u
1
u
2
u
3
u
4
=
0.4∆x
2
− (1 + 0.5c∆x)u
0
0.4∆x
2
0.4∆x
2
0.4∆x
2
− (1 + 4.5c∆x)u
5
(3)
By discretizing the differential equations we have created a set of linear algebraic equations that have
the symbolic form AU = b. To set up and solve the equations (3) in Matlab, we could type:
3
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.
Exercise 1: In this exercise, you will be using the above code to solve the rope BVP. You will also be
exhaustively checking that the code is correct.
(a) Copy the above code and paste it into a script m-file named exer1a.m. Execute exer1a to find
a solution of Equation (3). Please include the printed values of U as part of the lab summary.
(b) Verify that the values of U and b that you found satisfy at least one of the middle four equations in
(2). To do this, write a script m-file named exer1b.m and plug the values of c, ∆x, u
k
, k = 0, . . . , 5
into your chosen equation. Show the result is essentially zero.
Be careful! Lower-case c is upper-case C, ∆x is dx and u
k
, k = 0, 1, . . . 5 is U(1:6) in exer1a.
(c) Direct substitution into (3) shows that the function u(x) = 1 would be a solution if ρg = 0 and
u(5) = 1. 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.
(d) Direct substitution into (3) shows that the function u(x) = x would b e a solution if ρg = c,
u(1) = 0 and u(5) = 5. 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.
(e) Direct substitution into (3) shows that the function u(x) = x
2
would be a solution if ρg =
2 + 4cx, u(1) = 0) and u(5) = 25. 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 u(x) = x
2
and check that the
discrete solution is (exactly or to roundoff) correct.
(f) Now that you are confident that the code is correct, use exer1a.m to solve the unmodified BVP
(3). Plot U versus x. It should appear roughly parabolic, like a rope hanging from its ends, and
pass through (0, 1) and (5, 1.5) on its ends. Please include this plot with your summary.
4
剩余16页未读,继续阅读
资源评论
卷积神经网络
- 粉丝: 340
- 资源: 8460
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功