-
Notifications
You must be signed in to change notification settings - Fork 0
/
penta2.m
84 lines (73 loc) · 1.47 KB
/
penta2.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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
function [x]= penta2(y,a,b,c)
%
% Author: Kurt Annen [email protected]
% Date: 15/05/2004
% Internet: www.web-reg.de
%
% Solves the problem Ax=b when A is pentadiagonal and strongly nonsingular.
% This is much faster than x=A\y for large matrices.
%
% Reference: Späth, Helmuth "Numerik: Eine Einführung für Mathematiker und Informatiker"
% S. 110 . Vieweg-Verlag Braunschweig/Wiesbaden (1994)
%
% a = main diagonal
% b = 2. diagonal
% c = 3. diagonal
%
if nargin ~= 4
error('penta(A,y) requires four arguments');
end
[n] = length(a);
[m] = length(y);
[o] = length(b);
[p] = length(c);
if m ~= n
error('Inner matrix dimensions must agree');
end
c(m)=0;
c(m-1)=0;
b(m)=0;
if m ~= length(b) & m ~= length(c)
error('a,b,c must have the same dimension');
end
% a optimized algorithm
h1=0;
h2=0;
h3=0;
h4=0;
h5=0;
hh1=0;
hh2=0;
hh3=0;
hh4=0;
hh5=0;
z=0;
hb=0;
hc=0;
for i=1:m
z=a(i)-h4*h1-hh5*hh2;
hb=b(i);
hh1=h1;
h1=(hb-h4*h2)/z;
b(i)=h1;
hc=c(i);
hh2=h2;
h2=hc/z;
c(i)=h2;
a(i)=(y(i)-hh3*hh5-h3*h4)/z;
hh3=h3;
h3=a(i);
h4=hb-h5*hh1;
hh5=h5;
h5=hc;
end
h2=0;
h1=a(m);
y(m)=h1;
for i=m:-1:1
y(i)=a(i)-b(i)*h1-c(i)*h2;
h2=h1;
h1=y(i);
end
x=y;
return