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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
// Computation of residues
// 4.5
// Numerator and denominator coefficients
// are passed in decreasing powers of z(say)
function [res,pol,q] = respol(num,den)
len = length(num);
if num(len) == 0
num = num(1:len-1);
end
[resi,q] = pfe(num,den);
res = resi(:,2);
res = int(res) + (clean(res - int(res),1.d-04));
pol = resi(:,1);
pol = int(pol) + (clean(pol - int(pol),1.d-04));
endfunction;
///////////////////////////////////////////////////
// Partial fraction expansion
function [resid1,q] = pfe(num,den)
x = poly(0,'x');
s = %s;
exec flip.sci;
num2 = flip(num);
den2 = flip(den);
num = poly(num2,'s','coeff');
den = poly(den2,'s','coeff');
[fac,g] = factors(den);
polf = polfact(den);
n = 1;
r = clean(real(roots(den)),1.d-5);
//disp(r);
l = length(r);
r = gsort(r,'g','i');
r = [r; 0];
j = 1;
t1 = 1; q = [];
rr = 0;
rr1 = [0 0];
m = 1;
for i = j:l
if abs(r(i)- r(i+1)) < 0.01 then
r(i);
r(i+1);
n = n+1;
m = n;
//pause
t1 = i;
//disp('Repeated roots')
else
m = n;
//pause
n = 1;
end
i;
if n == 1 then
rr1 = [rr1; r(i) m];
//pause
end;
j = t1 + 1;
end;
rr2 = rr1(2:$,:);
[r1,c1] = size(rr2);
den1 = 1;
for i = 1:r1
den1 = den1 * ((s-rr2(i,1))^(rr2(i,2)));
end;
[rem,quo] = pdiv(num,den);
q = quo;
if quo ~= 0
num = rem;
end
tf = num/den1;
res1 = 0;
res3 = [s 0];
res5 = [0 0];
for i = 1:r1
j = rr2(i,2) + 1;
tf1 = tf; //strictly proper
k = rr2(i,2);
tf2 = ((s-rr2(i,1))^k)*tf1;
rr2(i,1);
res1 = horner(tf2,rr2(i,1));
res2 = [(s - rr2(i,1))^(rr2(i,2)) res1];
res4 = [rr2(i,1) res1];
res3 = [res3; res2];
res5 = [res5; res4];
res = res1;
for m = 2:j-1
k;
rr2(i,1);
tf1 = derivat(tf2)/factorial(m-1); //ith derivative
res = horner(tf1,rr2(i,1));
res2 = [(s - rr2(i,1))^(j-m) res];
res4 = [rr2(i,1) res];
res5 = [res5; res4];
res3 = [res3; res2];
tf2 = tf1;
end;
end;
resid = res3(2:$,:); //with s terms
resid1 = res5(2:$,:); //only poles(in decreasing no. of repetitions)
endfunction;
////////////////////////////////////////////////////////////
|
