-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcheckjac.m
121 lines (102 loc) · 3.07 KB
/
checkjac.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
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
113
114
115
116
117
118
119
120
121
numjacs
fm_call('i')
disp(' ')
gerr = max(abs(DAE.g));
ferr = max(abs(DAE.f));
disp(['g max err = ',num2str(gerr)])
disp(['f max err = ',num2str(ferr)])
if ferr > Settings.dyntol
disp('The following elements of the f vector are suspiciously high:')
disp(' ')
v = abs(DAE.f);
i = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h));
disp(['* ', Varname.uvars{i(h)}, ' -> ', num2str(u)])
end
disp(' ')
end
if gerr > Settings.dyntol
disp('The following elements of the g vector are suspiciously high:')
disp(' ')
v = abs(DAE.g);
i = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h));
disp(['* ', Varname.uvars{i(h) + DAE.n}, ' -> ', num2str(u)])
end
disp(' ')
end
disp(' ')
disp(['Fx abs err = ',num2str(max(max(abs(DAE.Fx-Fx))))])
disp(['Fy abs err = ',num2str(max(max(abs(DAE.Fy-Fy))))])
disp(['Gx abs err = ',num2str(max(max(abs(DAE.Gx-Gx))))])
disp(['Gy abs err = ',num2str(max(max(abs(DAE.Gy-Gy))))])
disp(' ')
if DAE.n
[i,j] = find(abs(DAE.Fx) <= Settings.dyntol);
Fxtmp = DAE.Fx + sparse(i,j,1,DAE.n,DAE.n);
[i,j] = find(abs(DAE.Fy) <= Settings.dyntol);
Fytmp = DAE.Fy + sparse(i,j,1,DAE.n,DAE.m);
[i,j] = find(abs(DAE.Gx) <= Settings.dyntol);
Gxtmp = DAE.Gx + sparse(i,j,1,DAE.m,DAE.n);
else
Fxtmp = 1;
Fytmp = ones(1,DAE.m);
Gxtmp = ones(DAE.m,1);
end
[i,j] = find(abs(DAE.Gy) <= Settings.dyntol);
Gytmp = DAE.Gy + sparse(i,j,1,DAE.m,DAE.m);
Fxerr = max(max(abs((DAE.Fx-Fx)./Fxtmp)));
Fyerr = max(max(abs((DAE.Fy-Fy)./Fytmp)));
Gxerr = max(max(abs((DAE.Gx-Gx)./Gxtmp)));
Gyerr = max(max(abs((DAE.Gy-Gy)./Gytmp)));
disp(['Fx rel err = ',num2str(Fxerr)])
disp(['Fy rel err = ',num2str(Fyerr)])
disp(['Gx rel err = ',num2str(Gxerr)])
disp(['Gy rel err = ',num2str(Gyerr)])
disp(' ')
if Fxerr > Settings.dyntol
disp('The following elements of the Fx matrix are suspiciously high:')
disp(' ')
v = abs((DAE.Fx-Fx)./Fxtmp);
[i, j] = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h), j(h));
disp(['* ', Varname.uvars{i(h)}, ' - ', Varname.uvars{j(h)}, ' -> ', num2str(u)])
end
disp(' ')
end
if Fyerr > Settings.dyntol
disp('The following elements of the Fy matrix are suspiciously high:')
disp(' ')
v = abs((DAE.Fy-Fy)./Fytmp);
[i, j] = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h), j(h));
disp(['* ', Varname.uvars{i(h)}, ' - ', Varname.uvars{j(h)+DAE.n}, ' -> ', num2str(u)])
end
disp(' ')
end
if Gxerr > Settings.dyntol
disp('The following elements of the Gx matrix are suspiciously high:')
disp(' ')
v = abs((DAE.Gx-Gx)./Gxtmp);
[i, j] = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h), j(h));
disp(['* ', Varname.uvars{i(h)+DAE.n}, ' - ', Varname.uvars{j(h)}, ' -> ', num2str(u)])
end
disp(' ')
end
if Gyerr > Settings.dyntol
disp('The following elements of the Gy matrix are suspiciously high:')
disp(' ')
v = abs((DAE.Gy-Gy)./Gytmp);
[i, j] = find(v > Settings.dyntol);
for h = 1:length(i)
u = v(i(h), j(h));
disp(['* ', Varname.uvars{i(h)+DAE.n}, ' - ', Varname.uvars{j(h)+DAE.n}, ' -> ', num2str(u)])
end
disp(' ')
end