Solving ODEs with unknown parameter using bvp5c (Updated)

Document technical information

Format pdf
Size 244.8 kB
First found Nov 13, 2015

Document content analysis

Category Also themed
Language
English
Type
not defined
Concepts
no text concepts found

Transcript

1
2
%% C l e a r workspace
clc ; clear all ; close all ;
3
Solving ODEs with unknown parameter using bvp5c (Updated)
4
5
6
7
Jaehyuck Park
[email protected]
8
9
%% C o n s t a n t s
r = 0.08;
Sigma = 0 . 1 5 ;
phi = 0 . 1 5 ;
C = 0.1;
N = 100;
10
11
April 26, 2014
12
13
%% Numerical s o l u t i o n
d f = @( t , x , X B ) [ x ( 2 ) ; 2 / ( 1 ˆ 2 ⇤ ( t+X B ( 1 ) ) ˆ 2 ) ⇤ ( r ⇤x ( 1 ) r ⇤ ( t+X B ( 1 ) ) ⇤x ( 2 ) (1 p h i ) ⇤ ( t+X B ( 1 ) C) ) ] ;
bc = @( ya , yb , X B ) [ ya ( 1 ) ; ya ( 2 ) ; yb ( 2 ) 1];
14
15
Problem
16
17
18
Using MATLAB bvp5c, solve the following ODE numerically and find XB that satisfies (2).
(1
)(X
1
C) +
2
19
20
2
2
00
0
X F (X) + rXF (X)
rF (X) = 0
(1)
21
22
where F (XB ) = 0 and F 0 (XB ) = 0.
% Mesh
x1 = [ 0 ; 0 ] ;
x i n t = l i n s p a c e ( 0 , 1 0 ,N) ;
s i n i t=b v p i n i t ( x i n t , x1 , 0 . 1 ) ;
(2)
23
% Solver
s o l=bvp5c ( df , bc , s i n i t ) ;
Sxint = deval ( sol , xint ) ;
X BT = s o l . p a r a m e t e r s ( 1 ) % show undermined p a r a m e t e r s
24
The first condition of (2) is so-called ’value-matching condition’.
25
26
27
Approach 1
% Plot
plot ( sol . x , sol . y (1 ,:) ) ;
t i t l e ( ’ S o l u t i o n t o ODE with unknown parameter X B ’ ) ; x l a b e l ( ’X ’ ) ; y l a b e l ( ’F(X) ’ )
Results:
1
The bvp5c function solves ODEs on the interval [a,b] subject to two-point boundary value conditions
X BT =
2
bc(y(a), y(b)) = 0
3
23.5057
The two points (starting and end points) for the boundary conditions are first defined to make mash for the
numerical ODE solver. So we need to type in bcfun:
sol = bvp5c(odefun, bcfun, solinit)
res = bcfun(ya, yb, parameters)
where ya and yb are column vectors corresponding to y(a) and y(b). However in this case, these points are not
given. Here is the trick. We make use of the following change of variable technique:
Y ⌘X
XB
G(Y ) = G(X
XB ) = F (X)
Then the equations (1) and (2) changes to
(1
)(Y + XB
C) +
1
2
2
(Y + XB )2 G00 (Y ) + r(Y + XB )G0 (Y )
rG(Y ) = 0
(3)
where G(0) = 0 and G0 (0) = 0.
(4)
Plus,
lim G0 (Y ) = 1
Y !1
(5)
Now we can solve the ODE using bvp5c letting XB as a parameter to be determined. Below shows the MATLAB
code. Here I assumed the constant values r = 0.08, = 0.15, = 0.15 and C = 0.1.
1
2
Approach 2
1
2
Also we can solve this Euler-Lagrange equation by hand. First we find homogeneous solutions y1 and y2 . If we
let yh = x↵ and substitute this to equation (1) we get
⇢
y1 = x
2
y2 = x r/a = x 2r/
Then using the MATLAB dsolve function, we find the particular solution yp .
✓
◆
1
C1 2
x ln x + C +
yp = 1 2
r
2
+
r
2
3
4
5
%% Symbolic s o l u t i o n
syms x a b C r a l p h a Sigma p h i A1 A2
6
7
8
9
10
11
12
13
14
15
Adding these three solutions, we get the general solution yg
%% C l e a r t h e workspace
clc ; clear all ; close all ;
16
ans = d s o l v e ( ’ (1 p h i ) ⇤ ( x C) +1/2⇤Sigma ˆ2⇤ x ˆ2⇤D2y+r ⇤x⇤Dy r ⇤y=0 ’ , x )
% Homogeneous s o l u t i o n s
y1 = x ;
y2 = xˆ( r / ( 1 / 2 ⇤ Sigma ˆ 2 ) ) ;
% Particular solution
yp = ( phi 1) / ( 1 / 2 ⇤ Sigmaˆ2+ r ) ⇤ (C+x⇤ l o g ( x )+C/ r ⇤1/2⇤ Sigma ˆ 2 )
% General s o l u t i o n
yg = A1⇤ y1+A2⇤ y2+yp ;
17
yg = A1 y1 + A2 y2 + yp
= A1 x + A2 x
2r/
2
+
1
2
1
2+r
✓
x ln x + C +
C1
r 2
2
18
◆
19
20
21
22
% Check i f yp and yg s a t i s f y t h e g i v e n ODE
r e s u l t p = (1 p h i ) ⇤ ( x C) +1/2⇤Sigma ˆ2⇤ x ˆ2⇤ d i f f ( yp , x , 2 )+r ⇤x⇤ d i f f ( yp , x , 1 ) r ⇤yp ;
simplify ( result p )
r e s u l t g = (1 p h i ) ⇤ ( x C) +1/2⇤Sigma ˆ2⇤ x ˆ2⇤ d i f f ( yg , x , 2 )+r ⇤x⇤ d i f f ( yg , x , 1 ) r ⇤ yg ;
simplify ( result g )
where A1 and A2 are constants to be determined by the boundary conditions.
Results:
Assumption 1. If we use boundary conditions
0
(1
then we get A1 = 2.08, A2 =
F (XB ) = 0 and F (XB ) = 0
1 2 2 00
)(XB C) +
XB F (XB ) = 0
2
(6)
(7)
1.02, and XB = 1.00.
1
ans =
2
3
C2⇤x + C3⇤ exp ( (2⇤ r ⇤ l o g ( x ) ) / Sigma ˆ 2 ) + ( 2 ⇤ (C + x⇤ l o g ( x ) ) ⇤ ( p h i
1 ) ) / ( Sigma ˆ2 + 2⇤ r ) + ( Sigma ˆ2⇤
x ˆ ( ( 2 ⇤ r ) / Sigma ˆ 2 ) ⇤ exp ( (2⇤ r ⇤ l o g ( x ) ) / Sigma ˆ 2 ) ⇤ ( p h i
1 ) ⇤ (C⇤ Sigma ˆ2 + 2⇤C⇤ r
2⇤ r ⇤x ) ) / ( r ⇤ (
Sigma ˆ2 + 2⇤ r ) ˆ 2 )
4
5
6
8
F (XB ) = 0 and F 0 (XB ) = 0
(8)
lim F 0 (X) = 1
(9)
X!1
yp =
7
Assumption 2. However, if we use boundary conditions
then there is no constant values A1 , A2 , and XB satisfying these since there is a term (x ln x) which has an
divergent derivative. Attached codes shows how to find the general symbolic solution.
( ( phi
9
1 ) ⇤ ( (C⇤ Sigma ˆ 2 ) / ( 2 ⇤ r ) + C + x⇤ l o g ( x ) ) ) / ( Sigma ˆ2/2 + r )
10
11
ans =
12
13
0
14
15
16
ans =
17
18
3
0
4

Similar documents

×

Report this document