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
|
(kill(all), load(simplex), 0);
0;
linear_program(matrix([1,1,-1,0],[2,-3,0,-1],[4,-5,0,0]), [1,1,6], [1,-2,0,0]);
[[13/2, 4, 19/2, 0], -3/2];
minimize_lp(x, [y>=x-1, y>=-x-1, y<=x+1, y<=1-x, y=x/2]);
[-(2/3), [x=-(2/3), y=-(1/3)]];
maximize_lp(x+y, [y<=-x/2+3, y<=-x+4], [x, y]);
[4, [x = 2, y = 2]];
maximize_lp(x+y, [y<=-x/2+3, y<=-x+4], all);
[4, [x = 2, y = 2]];
maximize_lp(x+y, [y<=-x/2+3, y<=-x+4]), nonnegative_lp=true;
[4, [x = 2, y = 2]];
maximize_lp(10*x1 - 57*x2 - 9*x3 - 24*x4,
[1/2*x1 - 11/2*x2 - 5/2*x3 + 9*x4 < 0,
1/2*x1 - 3+2*x2 - 1/2*x3 + x4 < 0,
x1 < 1],
all)$
[41/5,[x4=0,x3=1/5,x2=0,x1=1]];
minimize_lp(x, [x<1, y>0]);
"Problem not bounded!";
minimize_lp(x, [x<1, y>0], all);
[0,[y=0,x=0]];
minimize_lp(x, [x<1, y<-1], all);
"Problem not feasible!";
block([A,b,c],
A : read_matrix(file_search("Tests/afiro_A.csv"), 'csv),
b : read_list(file_search("Tests/afiro_b.csv"), 'csv),
c : read_list(file_search("Tests/afiro_c.csv"), 'csv),
is (abs(second(linear_program(A, b, c) + 464.7531428571429))<10^-4)
);
true;
block([A,b,c],
A : read_matrix(file_search("Tests/sc50a_A.csv"), 'csv),
b : read_list(file_search("Tests/sc50a_b.csv"), 'csv),
c : read_list(file_search("Tests/sc50a_c.csv"), 'csv),
is (abs(second(linear_program(A, b, c) + 64.5750770585645))<10^-4)
);
true;
minimize_lp(-x, [1e-9*x + y <= 1], [x,y])$
"Problem not bounded!"$
is(abs(first(minimize_lp(-x, [1e-9*x + y <= 1], [x,y])) + 1e9) < 1e-5), epsilon_lp=0 $
true;
minimize_lp(-x, [10^-9*x + y <= 1], [x,y])$
[- 1000000000, [y = 0, x = 1000000000]]$
/* Thanks to Michael Soegtrop
* https://sourceforge.net/p/maxima/mailman/message/36414030/
*
* Do max-norm interpolation of a polynomial (x^3) on an interval
* ([3/4,1]) with a quintic and sextic at 10 points. Using rational
* arithmetic, we should get x^3; using floating point arithmetic and the
* default value of epsilon_lp, we do for the quintic (1.0*x^3) but not
* the sextic.
*
*/
block([solutionpoly_lmad,f,sol5,sol6,sol5f,sol6f,checkerr],
local(solutionpoly_lmad,f,checkerr),
solutionpoly_lmad(f,degree,npoints,xmin,xmax,float_coerce):=block([xp,gep,gem,ineq,sol,yp,ydiff],
xp:float_coerce(makelist(xmin+(xmax-xmin)*(i-1)/(npoints-1),i,1,npoints)),
yp:makelist(apply(f,[xp[i]]),i,1,npoints),
ydiff:makelist(yp[i]-sum(a[k]*xp[i]^k,k,0,degree),i,1,npoints),
gep:makelist(umax>=ydiff[i],i,1,npoints),
gem:makelist(umax>=-ydiff[i],i,1,npoints),
ineq:append(gep,gem),
sol:second(minimize_lp(umax,ineq)),
/*[at(sum(a[i]*'x^i,i,0,degree),sol),at(umax,sol),sol]*/
at(umax,sol)
),
f(x):=x^3,
sol5:solutionpoly_lmad(f,5,10,3/4,1,identity),
sol6:solutionpoly_lmad(f,6,10,3/4,1,identity),
sol5f:solutionpoly_lmad(f,5,10,3/4,1,'float),
sol6f:solutionpoly_lmad(f,6,10,3/4,1,'float),
checkerr(u,v) := is(abs(u)<=v),
map(checkerr, [sol5,sol6,sol5f,sol6f,sol6f,sol6f],[0,0,1e-16,1e-7,1e-10,sol5f])
);
[true,true,true,true,false,false] $
(assume(c>4/3), declare(c,constant),
ratsimp(maximize_lp(x+y, [y<=-x/c+3, y<=-x+4], [x, y]) - [4,[x = c/(c-1),y = (3*c-4)/(c-1)]])), epsilon_lp=0 $
[0,[0 = 0,0 = 0]] $
(forget(c>4/3), remove(c,constant), kill(nonnegative_lp,return_input_lp), 0)$
0 $
|
