Lecture 30: Systems of Ordinary Differential

1
Lecture 30 Systems of Ordinary
Differential Equations and Runge-Kutta methods
Eulers method for systems
euler2.m
Program 6.2 Vector version of Euler method
Inputs interval a,b, initial vector y0, step
size h Output time steps t, solution y
Example usage t yeuler2(0 1,0
1,0.1) function t yeuler2(int,y0,h) t(1)int(
1) y(1,)y0 nround((int(2)-int(1))/h) for
i1n t(i1)t(i)h y(i1,)eulerstep(t(i),y
(i,),h) end function yeulerstep(t,y,h) one
step of the Euler method Input current time t,
current vector y, step size h Output the
approximate solution vector at time
th yyhydot(t,y) function zydot(t,y) z(1)
y(2)2-2y(1) z(2) y(1)-y(2)-ty(2)2
2
euler2test.m - test of function euler2. Runing
require file euler2.m Example usage
yeuler2(0 1,0 1,0.1) y1exactinline('t.ex
p(-2t)') y2exactinline('exp(-t)') t
yeuler2(0 1, 0 1, 0.01) plot(t,y(,1),t,y(
,2)) legend('y1','y2') y1exactvecy1exact(t)
y2exactvecy2exact(t) ii110 hvec0.1./2.ii
for i1length(hvec) t yeuler2(0 1,
0 1, hvec(i)) disp('h',num2str(hvec(i))
,' Error in eulers method for y(1)
', num2str(y(length(y))-y1exactvec(length(y1exact
vec)))) end
3
gtgt euler2test h0.05 Error in eulers method for
y(1) 0.0057103 h0.025 Error in eulers method
for y(1) 0.0028366 h0.0125 Error in eulers
method for y(1) 0.0014133 h0.00625 Error in
eulers method for y(1) 0.00070534 h0.003125
Error in eulers method for y(1)
0.00035234 h0.0015625 Error in eulers method
for y(1) 0.00017609 h0.00078125 Error in eulers
method for y(1) 8.8023e-005 h0.00039063 Error
in eulers method for y(1) 4.4006e-005 h0.0001953
1 Error in eulers method for y(1)
2.2002e-005 h9.7656e-005 Error in eulers method
for y(1) 1.1001e-005
4
Eulers method for system of 4 equations Orbital
motion around sun
orbit.m
5
Program 6.4 Plotting program for one-body
problem Inputs inta b time interval,
initial conditions ic x0 vx0 y0 vy0, x
position, x velocity, y pos, y vel, h
stepsize, p steps per point plotted Calls a
one-step method such as trapstep.m Example
usage orbit(0 100,0 1 2 0,0.01,5) function
zorbit(int,ic,h,p) nround((int(2)-int(1))/(ph))
plot n points x0ic(1)vx0ic(2)y0ic
(3)vy0ic(4) grab initial conds y(1,)x0
vx0 y0 vy0t(1)int(1) build y
vector set(gca,'XLim',-5 5,'YLim',-5
5,'XTick',-5 0 5,'YTick',... -5 0
5,'Drawmode','fast','Visible','on','NextPlot','ad
d') cla sunline('color','y','Marker','.','marke
rsize',25,... 'xdata',0,'ydata',0) drawnow hea
dline('color','r','Marker','.','markersize',25,..
. 'erase','xor','xdata',,'ydata',) taillin
e('color','b','LineStyle','-','erase','none',...
'xdata',,'ydata',) px,py,buttonginput(1)
include these three
lines px1,py1,buttonginput(1) to
enable mouse support y(1,)px px1-px py
py1-py 2 clicks set direction for k1n
for i1p t(i1)t(i)h
y(i1,)eulerstep(t(i),y(i,),h) end
y(1,)y(p1,)t(1)t(p1) set(head,'xdata',y(
1,1),'ydata',y(1,3)) set(tail,'xdata',y(2p,1),'
ydata',y(2p,3)) drawnow End
6
function yeulerstep(t,x,h) one step of the
Euler method yxhydot(t,x) function z
ydot(t,x) m23g1mg2m2gpx20py20 px1x(1)
py1x(3)vx1x(2)vy1x(4) distsqrt((px2-px1)2
(py2-py1)2) zzeros(1,4) z(1)vx1 z(2)(mg2(p
x2-px1))/(dist3) z(3)vy1 z(4)(mg2(py2-py1))/
(dist3)
7
gtgt orbit(0 100, 0 1 2 0,0.01,5)
Show animation
8
Orbital motion around sun but now with trapezoid
method
orbittrap.m
obrittrap. - using trapezoid method for
integration Inputs inta b time interval,
initial conditions ic x0 vx0 y0 vy0, x
position, x velocity, y pos, y vel, h
stepsize, p steps per point plotted Calls a
one-step method such as trapstep.m Example
usage orbittrap(0 100,0 1 2
0,0.01,5) function zorbit(int,ic,h,p) nround((i
nt(2)-int(1))/(ph)) plot n points
x0ic(1)vx0ic(2)y0ic(3)vy0ic(4) grab
initial conds y(1,)x0 vx0 y0 vy0t(1)int(1)
build y vector set(gca,'XLim',-5
5,'YLim',-5 5,'XTick',-5 0 5,'YTick',...
-5 0 5,'Drawmode','fast','Visible','on','NextPlo
t','add') cla sunline('color','y','Marker','.',
'markersize',25,... 'xdata',0,'ydata',0) drawno
w headline('color','r','Marker','.','markersize'
,25,... 'erase','xor','xdata',,'ydata',) ta
illine('color','b','LineStyle','-','erase','none'
,... 'xdata',,'ydata',) px,py,buttongin
put(1) include these three
lines px1,py1,buttonginput(1) to
enable mouse support y(1,)px px1-px py
py1-py 2 clicks set direction
9
for k1n for i1p t(i1)t(i)h
y(i1,)trapstep(t(i),y(i,),h) end
y(1,)y(p1,)t(1)t(p1) set(head,'xdata',y(
1,1),'ydata',y(1,3)) set(tail,'xdata',y(2p,1),'
ydata',y(2p,3)) drawnow end function y
trapstep(t,x,h) one step of the Trapezoid
Method z1ydot(t,x) gxhz1 z2ydot(th,g) yx
h(z1z2)/2 function z ydot(t,x) m23g1mg
2m2gpx20py20 px1x(1)py1x(3)vx1x(2)vy1
x(4) distsqrt((px2-px1)2(py2-py1)2) zzeros
(1,4) z(1)vx1 z(2)(mg2(px2-px1))/(dist3) z(
3)vy1 z(4)(mg2(py2-py1))/(dist3)
10
gtgt orbittrap(0 100, 0 1 2 0,0.01,5)
Show animation
11
Hodgkin-Huxley equations as example or 4th order
Runge-Kutta method use
12
Program 6.5 Hodgkin-Huxley equations Inputs
a b time interval, ic initial voltage v
and 3 gating variables, step size h Output
solution y Calls a one-step method such as
rk4step.m Example usage yhh(0 100,-65 0
0.3 0.6,0.05) function yhh(inter,ic,h) global
pa pb pulse inpinput('pulse start, end, muamps
in , e.g. 50 51 7 ') painp(1)pbinp(2)pu
lseinp(3) ainter(1)binter(2)nceil((b-a)/h)
plot n points in total y(1,)ic
initial conditions t(1)a for
i1n t(i1)t(i)h y(i1,)rk4step(t(i),y(i
,),h) end subplot(3,1,1) plot(a pa pa pb pb
b,0 0 pulse pulse 0 0) gridaxis(0 100 0
2pulse) ylabel('input pulse') subplot(3,1,2) pl
ot(t,y(,1))gridaxis(0 100 -100
100) ylabel('voltage (mV)') subplot(3,1,3) plot(
t,y(,2),t,y(,3),t,y(,4))gridaxis(0 100 0
1) ylabel('gating variables') legend('m','n','h')
xlabel('time (msec)')
hh.m
13
function yrk4step(t,w,h) one step of the
Runge-Kutta order 4 method s1ydot(t,w) s2ydot(t
h/2,whs1/2) s3ydot(th/2,whs2/2) s4ydot(t
h,whs3) ywh(s12s22s3s4)/6 function
z ydot(t,w) global pa pb pulse c1g1120g236
g30.3T(papb)/2lenpb-pa e0-65e150e2-77
e3-54.4 inpulse(1-sign(abs(t-T)-len/2))/2
square pulse input on interval pa,pb of pulse
muamps vw(1)mw(2)nw(3)hw(4) zzeros(1,4)
z(1)(in-g1mmmh(v-e1)-g2nnnn(v-e2)-g3(v
-e3))/c v v-e0 z(2)(1-m)(2.5-0.1v)/(exp(2.
5-0.1v)-1)-m4exp(-v/18) z(3)(1-n)(0.1-0.01v
)/(exp(1-0.1v)-1)-n0.125exp(-v/80) z(4)(1-h)
0.07exp(-v/20) - h/(exp(3-0.1v)1)
14
Current above threshold 6.9 50 51 7
gtgt hh(0 100,-65 0 0.3 0.6,0.05) pulse start,
end, muamps in , e.g. 50 51 7 50 51 7
15
Current below threshold 6.9 50 51 5
gtgt hh(0 100,-65 0 0.3 0.6,0.05) pulse start,
end, muamps in , e.g. 50 51 7 50 51 5