-
Notifications
You must be signed in to change notification settings - Fork 0
/
p25.m
51 lines (51 loc) · 1.58 KB
/
p25.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
% p25.m - stability regions for ODE formulas
% Adams-Bashforth:
clf, subplot('position',[.1 .56 .38 .38])
plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
z = exp(1i*pi*(0:200)/100); r = z-1;
s = 1; plot(r./s)
s = (3-1./z)/2; plot(r./s)
s = (23-16./z+5./z.^2)/12; plot(r./s)
axis([-2.5 .5 -1.5 1.5]), axis square, grid on
title Adams-Bashforth
% order 1
% order 2
% order 3
% Adams-Moulton:
subplot('position',[.5 .56 .38 .38])
plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
s = (5*z+8-1./z)/12; plot(r./s)
% order 3
s = (9*z+19-5./z+1./z.^2)/24; plot(r./s)
% order 4
s = (251*z+646-264./z+106./z.^2-19./z.^3)/720; plot(r./s)
% 5
d = 1-1./z;
s = 1-d/2-d.^2/12-d.^3/24-19*d.^4/720-3*d.^5/160; plot(d./s) % 6
axis([-7 1 -4 4]), axis square, grid on, title Adams-Moulton
% Backward differentiation:
subplot('position',[.1 .04 .38 .38])
plot([-40 40],[0 0]), hold on, plot([0 0],[-40 40])
r = 0; for i = 1:6, r = r+(d.^i)/i; plot(r), end
% orders 1-6
axis([-15 35 -25 25]), axis square, grid on
title('backward differentiation')
% Runge-Kutta:
subplot('position',[.5 .04 .38 .38])
plot([-8 8],[0 0]), hold on, plot([0 0],[-8 8])
w = 0; W = w; for i = 2:length(z)
% order 1
w = w-(1+w-z(i)); W = [W; w]; end, plot(W)
w = 0; W = w; for i = 2:length(z)
% order 2
w = w-(1+w+.5*w^2-z(i)^2)/(1+w); W = [W; w];
end, plot(W)
w = 0; W = w; for i = 2:length(z)
% order 3
w = w-(1+w+.5*w^2+w^3/6-z(i)^3)/(1+w+w^2/2); W = [W; w];
end, plot(W)
w = 0; W = w; for i = 2:length(z)
% order 4
w = w-(1+w+.5*w^2+w^3/6+w.^4/24-z(i)^4)/(1+w+w^2/2+w.^3/6);
W = [W; w]; end, plot(W)
axis([-5 2 -3.5 3.5]), axis square, grid on, title Runge-Kutta