บทที่ 10
ปัญหาค่าขอบรอยค่อ (Boundary Value Problems)
10.1 ส่วนนำ
การติวเสริมนี้ออกแบบขึ้นมาเพื่อให้เราคุ้นเคยวิธีการยิงแก้ปัญหาค่าเงื่อนไขขอบ นี่เป็นองค์ประกอบหลักของคอมพิวเตอร์ของสุขภาพ
ปัญหา 10.1. แต่แรกสุดก่อนอื่นเราจะมองไปที่วิธีการแก้ปัญหาสมการทั่วไป
10.2 แก้สมการ
10.2.1 พล็อตกราฟ(Plotting)
สมมุติว่าเราต้องการหาค่าของ x = (x1, x2) กล่าวคือ
k(x1^2 + x2 ^2) -2 = 0
k
sin(x1x2) = 0 (10.1)
สำหรับบางค่าของ k.
บางครั้งสิ่งแรกที่เราควรทำคือพล็อตฟังก์ชันเหล่านี้ เราจำต้องกำหนดฟังก์ชันให้ ดังนั้นเราสามารถใช้รูทีนการพล็อตของไซแลบที่เป็นมาตรฐาน จึงแนะนำให้กำหนด 2 ฟังก์ชันสำหรับ 2 องค์ประกอบของฟังก์ชัน
function z = eqn1(x,y)
function z = eqn1(x,y)
z =
k*(x**2+y**2) - 2;
endfunction
function
z = eqn2(x,y)
z =
k*sin(x*y);
endfunction
เราสามารถพล็อตฟังก์ชันโดยใช้ fplot3d หรือฟังก์ชัน fcontour2d
ให้ทดลอง
k =
1;
x =
-4:0.1:4;
y =
x;
xbasc();
fplot3d(x,y,eqn2)
Notice that we have set the value of k and it is being passed through to the
function
eqn1 implicitly.
Contour
plots can sometimes be easier to work with. A contour plot for the
first
function can be obtained via
xbasc();
fcontour2d(x,y,eqn1,[0 5 10 20 25]);
Not
that we have asked for the contour lines for 0 up to 25 in steps of 5. Make
sure
you have a look at the plot.
Now
let us over plot with the contours of the second function.
fcontour2d(x,y,eqn2,
[0,1]);
The
overall plot is quite complicated (at least to me). Perhaps you might like to
plot
both functions separately to get an idea of the “shape” of both. But with
both
contour plots on the same graph, us the 2d zoom facility of Scilab to get
an
approximate solution (3 sig figs) of equation (10.1) (for k = 1).
10.2.2
Equation Solver
It
is useful to automate the solution of equations. We have studied the use of
Newton’e
method earlier in the course. Scilab provides a procedure to solve
general
equations. The procedure is called fsolve. This procedure requires an
initial
guess and the equations written as a Scilab function. For our problem
we
would produce a Scilab function as such
function
y = f(x)
y(1)
= eqn1(x(1),x(2));
y(2)
= eqn2(x(1),x(2));
endfunction
where
the functions eqn1 and eqn2 were defined in the last section. Of course
you
could define the function without using eqn1 and eqn2. Unfortunately the
fsolve
function expects a very particular form, which is different than the calling
form
for the plotting routines.
We
would store this in a file and use exec to load the function. We could
then
solve the equation using fsolve, as such
k=1;
[x,y]
= fsolve([0;0], f)
Computational
Science, Scilab Tutorials 69
Does
the solution provided by fsolve match the one you discovered manually?
It
might be useful to get some feedback about what fsolve is doing. One way
is
to provide some output from the function. Redefine the function f as such:
function
y = f(x)
y(1)
= eqn1(x(1),x(2));
y(2)
= eqn2(x(1),x(2));
printf(’x(1)=%15e
x(2)=%15e y(1)=%15e y(2)=%15e\n’,x(1),x(2),y(1),y(2))
endfunction
Take
note of the printf function. This is the formatted print command (like
the
corresponding C procedure). The formatting is accomplished using the codes
like
%15e which will produce a decimal number using 15 characters.
You
might even want to produce graphical output from your f function.
10.3
Two Point Boundary Problem
Consider
the two-point boundary value problem
y′′ = 10y3 + 3y + t2, 0 ≤ t ≤ 1,
with
boundary conditions
y(0)
= 0, y(1) = 1,
We
will try to use the shooting method to solve this problem.
To
use the shooting method we must first be able to solve the initial value
problem.
Then we must be able to test many initial conditions to find our required
value.
•
First setup a function which defines the differential equation using the call-
ing
sequence
function
udot = shootode(t,u)
•
Check your function by solving the ode with initial conditions y(0) = 0,
y′(0) = 1.
• Plot the solution. Mine looks like
the following
fshoot(1)
0 0.1 0.2 0.3
0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
1
2
3
4
5
6
• What is the
value of y at 1? (I got about 5.1).
• Take your ode
solver (with your plotting routines) and produce a function
with calling
sequence
function y =
fshoot(s)
which takes an
initial value for the slope at 0 and returns the value of
y(1) − 1.
For instance you should approximately get fshoot(1.0) should
return
approximately 4.1.
• Finally use
the fsolve function to find an approximation to the boundary
value problem.
Here is the
graphical result I obtained when running fsolve(1,fshoot).
Computational
Science, Scilab Tutorials 71
0 0.1 0.2 0.3
0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
1
2
3
4
5
6
10.4 Exercise
Solve the
two-point boundary value problem
y′′ = 8y3 + 3y − 10,
0 ≤ t ≤ 1,
with boundary
conditions
y(0) = 0, y(1)
= 2,
Hint: You might
have to hunt for a good starting value for the slope, some-
where between 2
and 3.9.
What is the
slope at 0 for the function which solve this boundary value prob-
lem?
ไม่มีความคิดเห็น:
แสดงความคิดเห็น