-
Notifications
You must be signed in to change notification settings - Fork 0
/
solve2dmonot.m
125 lines (115 loc) · 4.92 KB
/
solve2dmonot.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
122
123
124
125
# Принимает две функции F1(x, y) и F2(x, y),
# монотонно возрастающие по x и y,
# решает систему F1(x, y) = 0, F2(x, y) = 0,
# x0, y0 - начальное приближение
# x_eps - допустимая точность x, y
# F_eps - допустимая точность F
# R - максимальное значение |x| и |y|
# Возвращает матрицу N_roots x 2 из корней системы
# в порядке возрастания x
function sol = solve2dmonot (F1, F2, x0, y0, x_eps, F_eps, R)
sol = [];
# пускай приближение принадлежит линии {F1(x, y) = 0}
x00 = solve1dmonot(@(xx) F1(xx, y0), x0, x_eps, F_eps, R);
y00 = y0;
# вдруг мы попали в точности на корень
[x00, y00] = skip_root_right(F1, F2, x00, y00, 2 * x_eps, F_eps, R);
# пусть x возрастает, y убывает
x = x00;
y = y00;
while (abs(x) < R && abs(y) < R)
# идем дальше, чтобы перейти к поиску следующего корня
[x, y] = skip_root_right(F1, F2, x, y, 2 * x_eps, F_eps, R);
# приближение принадлежит линии {F1(x, y) = 0}
while (abs(x) < R && abs(y) < R && !is_root(F1, F2, x, y, x_eps, F_eps))
# F1(x, y) == 0
if (F2(x, y) > 0) # y будет убывать
y = solve1dmonot(@(yy) F2(x, yy), y, x_eps, F_eps, R);
# F2(x, y) == 0
else # x будет возрастать
x = solve1dmonot(@(xx) F2(xx, y), x, x_eps, F_eps, R);
# F2(x, y) == 0
endif
# F2(x, y) == 0
if (F1(x, y) > 0) # y будет убывать
y = solve1dmonot(@(yy) F1(x, yy), y, x_eps, F_eps, R);
# F1(x, y) = 0
else # x будет возрастать
x = solve1dmonot(@(xx) F1(xx, y), x, x_eps, F_eps, R);
# F1(x, y) == 0
endif
endwhile
if (is_root(F1, F2, x, y, x_eps, F_eps))
# добавляем найденную пару (x, y) к множеству корней
sol = [ sol ; x, y ];
endif
endwhile
# теперь пусть x убывает, y возрастает
x = x00;
y = y00;
while (abs(x) < R && abs(y) < R)
# идем дальше, чтобы перейти к поиску следующего корня
[x, y] = skip_root_left(F1, F2, x, y, 2 * x_eps, F_eps, R);
# приближение принадлежит линии {F1(x, y) = 0}
while (abs(x) < R && abs(y) < R && !is_root(F1, F2, x, y, x_eps, F_eps))
# F1(x, y) == 0
if (F2(x, y) < 0) # y будет возрастать
y = solve1dmonot(@(yy) F2(x, yy), y, x_eps, F_eps, R);
# F2(x, y) == 0
else # x будет убывать
x = solve1dmonot(@(xx) F2(xx, y), x, x_eps, F_eps, R);
# F2(x, y) == 0
endif
# F2(x, y) == 0
if (F1(x, y) < 0) # y будет возрастать
y = solve1dmonot(@(yy) F1(x, yy), y, x_eps, F_eps, R);
# F1(x, y) = 0
else # x будет убывать
x = solve1dmonot(@(xx) F1(xx, y), x, x_eps, F_eps, R);
# F1(x, y) == 0
endif
endwhile
if (is_root(F1, F2, x, y, x_eps, F_eps))
# добавляем найденную пару (x, y) к множеству корней
sol = [ x, y ; sol ];
endif
endwhile
endfunction
function [x, y] = skip_root_right (F1, F2, x0, y0, x_eps, F_eps, R)
x = x0;
y = y0;
while (is_root(F1, F2, x, y, x_eps, F_eps))
while (is_root(F1, F2, x, y, x_eps, F_eps))
x += x_eps;
y -= x_eps;
endwhile
# пускай приближение принадлежит линии {F1(x, y) = 0}
x = solve1dmonot(@(xx) F1(xx, y), x, x_eps, F_eps, R);
endwhile
# теперь is_root(F1, F2, x, y, x_eps, F_eps) == false
endfunction
function [x, y] = skip_root_left (F1, F2, x0, y0, x_eps, F_eps, R)
x = x0;
y = y0;
while (is_root(F1, F2, x, y, x_eps, F_eps))
while (is_root(F1, F2, x, y, x_eps, F_eps))
x -= x_eps;
y += x_eps;
endwhile
# пускай приближение принадлежит линии {F1(x, y) = 0}
x = solve1dmonot(@(xx) F1(xx, y), x, x_eps, F_eps, R);
endwhile
# теперь is_root(F1, F2, x, y, x_eps, F_eps) == false
endfunction
# проверяет, что x и y находятся вблизи кривых {F1 = 0} и {F2 = 0}
function ret = is_root (F1, F2, x, y, x_eps, F_eps)
F1_val = F1(x, y);
F2_val = F2(x, y);
s1 = sign(F1_val);
s2 = sign(F2_val);
is_root_F1 = (sign(F1(x - x_eps, y)) != s1) || (sign(F1(x + x_eps, y)) != s1) ...
|| (sign(F1(x, y - x_eps)) != s1) || (sign(F1(x, y + x_eps)) != s1);
is_root_F2 = (sign(F2(x - x_eps, y)) != s2) || (sign(F2(x + x_eps, y)) != s2) ...
|| (sign(F2(x, y - x_eps)) != s2) || (sign(F2(x, y + x_eps)) != s2);
ret = is_root_F1 && is_root_F2 || abs(F1_val) < F_eps && abs(F2_val) < F_eps;
endfunction