Vpasolve can not find a solution that I know it exist

조회 수: 5 (최근 30일)
clear
clc
%% Parameters and Assumptions
syms x11 x12 x21 x22 positive % Actions
syms a11 a12 a21 a22 positive % Marginal effects
syms c1 c2 positive % Costs
syms W1 W2 positive % Budget constraint
syms epsilon1 epsilon2 positive % Inital voters (can be seen as partisian voters)
syms b1 b2 positive % Total effects
syms V positive % Prize of winning
assume(a11 > 0 & a11 < 1)
assume(a12 > 0 & a12 < 1)
assume(a21 > 0 & a21 < 1)
assume(a22 > 0 & a22 < 1)
%% Theoretical Model
% Votes gain functions
f1 = b1*(x11^a11)*(x12^a12);
f2 = b2*(x21^a21)*(x22^a22);
F1 = f1 + epsilon1;
F2 = f2 + epsilon2;
% Payoff - Probability of winning.
P1 = (F1/(F1 + F2))*V;
P2 = (F2/(F1 + F2))*V;
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
isAlways(DP11 > 0)
ans = logical
1
isAlways(DP21 > 0)
ans = logical
1
% The probability of winning is strictly increasing in activity one (and also two).
% Therefore candidates spend all budget.
% Budget constraints
BC1 = W1 == x11*c1 + x12*c2;
BC2 = W2 == x21*c1 + x22*c2;
%% Theoretical Solution
BC1 = isolate(BC1,x12);
BC2 = isolate(BC2,x22);
P1 = subs(P1, lhs(BC1), rhs(BC1));
P2 = subs(P2, lhs(BC2), rhs(BC2));
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
eqns1 = [DP11 == 0];
vars1 = [x11];
S11 = solve(eqns1,vars1,'Real',true);
eqns2 = [DP21 == 0];
vars2 = [x21];
S21 = solve(eqns2,vars2,'Real',true);
% S11 and S21 are the equilibrium efforts of activity one for candidates one and two.
S12 = (W1 - S11*c1)/c2;
S22 = (W2 - S21*c1)/c2;
S12 = simplify(S12);
S22 = simplify(S22);
%% Equations
% First order conditions
Eq1 = ((a11*W1)/(c1*(a11+a12))) - x11 == 0;
Eq2 = ((a12*W2)/(c2*(a21+a22))) - x12 == 0;
% Election result (it is number of votes)
Eq3 = f1 - 1189313 == 0;
Eq4 = f2 - 1152271 == 0;
% Parameter restrictions
Eq5 = a11 + a12 - 1 == 0;
Eq6 = a21 + a22 - 1 == 0;
%% Substitution
% Candidate number one is LLP and candidate number two is DM.
Data = [28.752/c1,0.911/c2,29.664,31.012/c1,2.490/c2,33.502,1022013,1143681,1,1];
% The budget is in millions of UYU, epsilon1 and epsilon2 is number of inicial voters.
% Since do not have enough equations for parameters, assume b1 = b2 = 1.
% Last two digits of the vector "Data" reflect that.
E1 = subs(Eq1,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data);
E2 = subs(Eq2,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data);
E3 = subs(Eq3,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data);
E4 = subs(Eq4,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data);
E5 = a11 + a12 -1 == 0;
E6 = a21 + a22 -1 == 0;
%% Calibration
X0 = [0 1; 0 1; 0 1; 0 1; 0 Inf; 0 Inf];
S = vpasolve([E1,E2,E3,E4,E5,E6],[a11, a12, a21, a22, c1, c2],X0)
S = struct with fields:
a11: [0x1 sym] a12: [0x1 sym] a21: [0x1 sym] a22: [0x1 sym] c1: [0x1 sym] c2: [0x1 sym]
I know there is a solution to this system of equations but Matlab does not find it and I can see what is wrong on the code.

채택된 답변

Torsten
Torsten 2024년 2월 14일
편집: Torsten 2024년 2월 14일
Maybe it must read
Eq2 = ((a21*W2)/(c2*(a21+a22))) - x12 == 0;
instead of
Eq2 = ((a12*W2)/(c2*(a21+a22))) - x12 == 0;
?
clear
clc
%% Parameters and Assumptions
syms x11 x12 x21 x22 positive % Actions
syms a11 a12 a21 a22 positive % Marginal effects
syms c1 c2 positive % Costs
syms W1 W2 positive % Budget constraint
syms epsilon1 epsilon2 positive % Inital voters (can be seen as partisian voters)
syms b1 b2 positive % Total effects
syms V positive % Prize of winning
assume(a11 > 0 & a11 < 1)
assume(a12 > 0 & a12 < 1)
assume(a21 > 0 & a21 < 1)
assume(a22 > 0 & a22 < 1)
%% Theoretical Model
% Votes gain functions
f1 = b1*(x11^a11)*(x12^a12);
f2 = b2*(x21^a21)*(x22^a22);
F1 = f1 + epsilon1;
F2 = f2 + epsilon2;
% Payoff - Probability of winning.
P1 = (F1/(F1 + F2))*V;
P2 = (F2/(F1 + F2))*V;
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
% The probability of winning is strictly increasing in activity one (and also two).
% Therefore candidates spend all budget.
% Budget constraints
BC1 = W1 == x11*c1 + x12*c2;
BC2 = W2 == x21*c1 + x22*c2;
%% Theoretical Solution
BC1 = isolate(BC1,x12);
BC2 = isolate(BC2,x22);
P1 = subs(P1, lhs(BC1), rhs(BC1));
P2 = subs(P2, lhs(BC2), rhs(BC2));
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
eqns1 = [DP11 == 0];
vars1 = [x11];
S11 = solve(eqns1,vars1,'Real',true);
eqns2 = [DP21 == 0];
vars2 = [x21];
S21 = solve(eqns2,vars2,'Real',true);
% S11 and S21 are the equilibrium efforts of activity one for candidates one and two.
S12 = (W1 - S11*c1)/c2;
S22 = (W2 - S21*c1)/c2;
S12 = simplify(S12);
S22 = simplify(S22);
%% Equations
% First order conditions
Eq1 = ((a11*W1)/(c1*(a11+a12))) - x11 == 0;
Eq2 = ((a21*W2)/(c2*(a21+a22))) - x12 == 0;
% Election result (it is number of votes)
Eq3 = f1 - 1189313 == 0;
Eq4 = f2 - 1152271 == 0;
% Parameter restrictions
Eq5 = a11 + a12 - 1 == 0;
Eq6 = a21 + a22 - 1 == 0;
%% Substitution
% Candidate number one is LLP and candidate number two is DM.
Data = [28.752/c1,0.911/c2,29.664,31.012/c1,2.490/c2,33.502,1022013,1143681,1,1];
% The budget is in millions of UYU, epsilon1 and epsilon2 is number of inicial voters.
% Since do not have enough equations for parameters, assume b1 = b2 = 1.
% Last two digits of the vector "Data" reflect that.
E1 = subs(Eq1,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E1 = 
E2 = subs(Eq2,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E2 = 
E3 = subs(Eq3,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E3 = 
E4 = subs(Eq4,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E4 = 
E5 = a11 + a12 -1 == 0
E5 = 
E6 = a21 + a22 -1 == 0
E6 = 
syms c1s c2s
sola = solve([E1,E2,E5,E6],[a11 a12 a21 a22]);
E3 = subs(E3,[a11 a12 a21 a22],[sola.a11 sola.a12 sola.a21 sola.a22]);
E4 = subs(E4,[a11 a12 a21 a22],[sola.a11 sola.a12 sola.a21 sola.a22]);
E3 = E3 + 1189313;
E3 = simplify(log(lhs(E3))) == log(rhs(E3));
E3 = subs(E3,[log(c1) log(c2)],[c1s c2s]);
E4 = E4 + 1152271;
E4 = simplify(log(lhs(E4))) == log(rhs(E4));
E4 = subs(E4,[log(c1) log(c2)],[c1s c2s]);
solcs = solve([E3 E4],[c1s c2s]);
a11 = vpa(sola.a11)
a11 = 
0.96925566343042071197411003236246
a12 = vpa(sola.a12)
a12 = 
0.03074433656957928802588996763754
a21 = vpa(sola.a21)
a21 = 
0.027192406423497104650468628738583
a22 = vpa(sola.a22)
a22 = 
0.97280759357650289534953137126142
c1 = vpa(exp(solcs.c1s))
c1 = 
0.00002339002395546612021028972551146
c2 = vpa(exp(solcs.c2s))
c2 = 
0.0000021694431095803284252266058790349

추가 답변 (1개)

Walter Roberson
Walter Roberson 2024년 2월 14일
편집: Walter Roberson 2024년 2월 14일
There is no solution in the reals.
clear
clc
%% Parameters and Assumptions
syms x11 x12 x21 x22 positive % Actions
syms a11 a12 a21 a22 positive % Marginal effects
syms c1 c2 positive % Costs
syms W1 W2 positive % Budget constraint
syms epsilon1 epsilon2 positive % Inital voters (can be seen as partisian voters)
syms b1 b2 positive % Total effects
syms V positive % Prize of winning
assume(a11 > 0 & a11 < 1)
assume(a12 > 0 & a12 < 1)
assume(a21 > 0 & a21 < 1)
assume(a22 > 0 & a22 < 1)
%% Theoretical Model
% Votes gain functions
f1 = b1*(x11^a11)*(x12^a12);
f2 = b2*(x21^a21)*(x22^a22);
F1 = f1 + epsilon1;
F2 = f2 + epsilon2;
% Payoff - Probability of winning.
P1 = (F1/(F1 + F2))*V;
P2 = (F2/(F1 + F2))*V;
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
isAlways(DP11 > 0)
ans = logical
1
isAlways(DP21 > 0)
ans = logical
1
% The probability of winning is strictly increasing in activity one (and also two).
% Therefore candidates spend all budget.
% Budget constraints
BC1 = W1 == x11*c1 + x12*c2;
BC2 = W2 == x21*c1 + x22*c2;
%% Theoretical Solution
BC1 = isolate(BC1,x12);
BC2 = isolate(BC2,x22);
P1 = subs(P1, lhs(BC1), rhs(BC1));
P2 = subs(P2, lhs(BC2), rhs(BC2));
DP11 = diff(P1,x11);
DP21 = diff(P2,x21);
eqns1 = [DP11 == 0];
vars1 = [x11];
S11 = solve(eqns1,vars1,'Real',true);
eqns2 = [DP21 == 0];
vars2 = [x21];
S21 = solve(eqns2,vars2,'Real',true);
% S11 and S21 are the equilibrium efforts of activity one for candidates one and two.
S12 = (W1 - S11*c1)/c2;
S22 = (W2 - S21*c1)/c2;
S12 = simplify(S12);
S22 = simplify(S22);
%% Equations
% First order conditions
Eq1 = ((a11*W1)/(c1*(a11+a12))) - x11 == 0;
Eq2 = ((a12*W2)/(c2*(a21+a22))) - x12 == 0;
% Election result (it is number of votes)
Eq3 = f1 - 1189313 == 0;
Eq4 = f2 - 1152271 == 0;
% Parameter restrictions
Eq5 = a11 + a12 - 1 == 0;
Eq6 = a21 + a22 - 1 == 0;
%% Substitution
% Candidate number one is LLP and candidate number two is DM.
Data = [28.752/c1,0.911/c2,29.664,31.012/c1,2.490/c2,33.502,1022013,1143681,1,1];
% The budget is in millions of UYU, epsilon1 and epsilon2 is number of inicial voters.
% Since do not have enough equations for parameters, assume b1 = b2 = 1.
% Last two digits of the vector "Data" reflect that.
E1 = subs(Eq1,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E1 = 
E2 = subs(Eq2,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E2 = 
E3 = subs(Eq3,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E3 = 
E4 = subs(Eq4,[x11,x12,W1,x21,x22,W2,epsilon1,epsilon2,b1,b2],Data)
E4 = 
E5 = a11 + a12 -1 == 0
E5 = 
E6 = a21 + a22 -1 == 0
E6 = 
eqns = [E1, E2, E3, E4, E5, E6].';
partial_a11 = solve(eqns(1), a11, 'returnconditions', true)
partial_a11 = struct with fields:
a11: (599*a12)/19 parameters: [1×0 sym] conditions: a12 < 19/599
eqns2 = subs(eqns(2:end), a11, partial_a11.a11)
eqns2 = 
partial_a12 = solve(eqns2(4), a12)
partial_a12 = 
eqns3 = subs(eqns2([1 2 3 5]), a12, partial_a12)
eqns3 = 
partial_a21 = solve(eqns3(4), a21)
partial_a21 = 
eqns4 = subs(eqns3([1 2 3]), a21, partial_a21)
eqns4 = 
%partial_c2 = Inf; %by inspection
%at this point we have something divided by c2, raised to the power of a22.
%With c2 being inf to satisfy 3677/(30900*c2) == 0, normally
%(249/(100*c2))^a22 would give 0^a22 which would give 0, and 0 times something minus something cannot give 0.
% But if a22 is 0 then it would be 0^0 which would give 1, and we survive
% anther step
partial_a22 = 0; %by inspection
eqns5 = subs(eqns4(2:end), a22, partial_a22)
eqns5 = 
partial_c1 = solve(eqns5(2), c1)
partial_c1 = 
eqns6 = subs(eqns5(1), c1, partial_c1)
eqns6 = 
partial_c2 = solve(eqns6, c2)
partial_c2 = 
subs(eqns4(1), c2, partial_c2)
ans = 
  댓글 수: 1
Roque Ivan Acosta Bellizzi
Roque Ivan Acosta Bellizzi 2024년 2월 14일
Hi, you are right. Like @Torsten suggests, there is a mistake in the code. Thanks.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Equation Solving에 대해 자세히 알아보기

제품


릴리스

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by