i want to solve a 10 equation, 10 states with fsolve , but it does not give me the answer.

HI Guys. Actually I want t o test that is there a 4x4 matrix named P which is positive and A1'*P+P1*A1 and also A2'*P+P1*A2 are negative matrixes?

May you please guide me how can I test/ do that?

Firstly I try to solve a 2x2 matrix P with Fsolve in a way that A1'*P+P*A1=-eye(2,2) , and see is there any positive p that gives me this result, but actually the bellow code that I wrote does not give me any solution. Why?

And I will be thankful if any one can help me in my problem in 4x4 matrix A1,A2 that I write in my mfile.

clear all;
clc;
close all;
A1=[3.4*10^-15 7.6*10^-6;1 4];
%A1=[3.4*10^-15 7.6*10^-6 0 0;-4*10^-16 -6.5*10^-5 0 0;2.3*10^-16 0.00059 -0.1426 0;-8.159*10^-14 -0.05 18.21 0.08];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.977*10^-7 -1.008*10^-5 0 0;9.24*10^-8 0.00045 -0.127 0;-8.55*10^-5 -0.10 65.217 -0.0833];
syms p1
syms p2
syms p3
syms p4
syms p5
syms p6
syms p7
syms p8
syms p9
syms p10
%P=[p1 p2 p3 p4;p2 p5 p6 p7;p3 p6 p8 p9;p4 p7 p9 p10];
P=[p(1) p(2);p(2) p(3)];
p(1)>0;
det(P)>0;
function F = myfun(p)
F=A1'*P+P*A1+eye(2,2)=0;
p0=[1 1;1 1];
[p,fval]=solve(@myfun,p0)

[John D'Errico]

Note that your subject title says that you want to use fsolve, not solve. solve and fsolve are different things.

[Walter Roberson]

When you say a "positive" matrix, do you mean "positive definite", or "positive semi-definite", or just that all the values are greater than 0 ?

[Walter Roberson]

azam ghamari comments to John D'Errico

Ok, I want to solve firstly this problem (fsolve) and then my initial question. Thanks

[Walter Roberson]

azam ghamari comments to me:

yes, I mean positive definite.

Answer 1

[azam ghamari]

thank you, but no it gives me this error:

Error: File: Untitled2.m Line: 8 Column: 1

Function definitions are not permitted in this context.

and also after solving this, do you have any suggestion for the other part ?" Actually I want t o test that is there a 4x4 matrix named P which is positive and A1'*P+P1*A1 and also A2'*P+P1*A2 are negative matrixes.", I brout A1,A2 also in my mfile.

Thank you

[John D'Errico]

Please learn to use comments instead of adding multiple answers just ot make a response.

Answer 2

[Walter Roberson]

A1=[3.4*10^-15 7.6*10^-6;1 4];
%A1=[3.4*10^-15 7.6*10^-6 0 0;-4*10^-16 -6.5*10^-5 0 0;2.3*10^-16 0.00059 -0.1426 0;-8.159*10^-14 -0.05 18.21 0.08];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.977*10^-7 -1.008*10^-5 0 0;9.24*10^-8 0.00045 -0.127 0;-8.55*10^-5 -0.10 65.217 -0.0833];

P = sym('p', [1 3]);

sol = solve(myfun(P,A1));
p = vpa([sol.p1, sol.p2, sol.p3]);
disp(p)

function F = myfun(p, A1)
P=[p(1) p(2);p(2) p(3)];
F=A1'*P+P*A1+eye(2,2);
end

[Walter Roberson]

Before R2016b, you would need to have the myfun part in a separate myfun.m file.

[azam ghamari]

Actually I write this code and the myfun part in another mfile. but it gives me just this:

p1 =

p1

it doesn't solve!

clear all;
clc;
close all;
A1=[3.4*10^-15 7.6*10^-6;1 4];
%A1=[3.4*10^-15 7.6*10^-6 0 0;-4*10^-16 -6.5*10^-5 0 0;2.3*10^-16 0.00059 -0.1426 0;-8.159*10^-14 -0.05 18.21 0.08];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.977*10^-7 -1.008*10^-5 0 0;9.24*10^-8 0.00045 -0.127 0;-8.55*10^-5 -0.10 65.217 -0.0833];
syms p1
syms p2
syms p3
syms p4
syms p5
syms p6
syms p7
syms p8
syms p9
syms p10
syms P
%P=[p1 p2 p3 p4;p2 p5 p6 p7;p3 p6 p8 p9;p4 p7 p9 p10];
P=[p(1) p(2);p(2) p(3)];
p(1)>0;
det(P)>0;
[P,fval]=fsolve(@myfun,p0)

[Walter Roberson]

Your p0 needs to be quite close to a solution in order for the fsolve() to not fail. For example, [28000 -1/2 -1/10]

There is an exact solution for myfun with that A1, and the solution is unique. However, the P that results from it is not positive definite.

[Walter Roberson]

azam ghamari comments to me:

from which initial condition you get this result? I entered these, but i give me this error: Undefined function 'p' for input arguments of type 'double'. Error in p1 (line 19) P=[p(1) p(2);p(2) p(3)]; Actually with every initial condition, it gives me this error. Thanks

[Walter Roberson]

azam, do not use flags to comment on postings -- not unless the posting is abusive or needs to be brought to the attention of the people who run the site for some reason.

[Walter Roberson]

myfun([28000 -1/2 -1/10], A1)

The exact solution is

[ 213421363591200938830046067085490690614800743154049565589504/763295569657123614987662768294787948506729929111387655, -10904222444405555534397559487062678279629900156502016/21808444847346388999647507665565369957335140831753933, -21808279103165234035220121295486902434609061492848333/174467558778771111997180061324522959658681126654031464]

which is approximately

[279605.138658241 -0.500000000950657 -0.124999049999998]

[azam ghamari]

Walter, myfun([28000 -1/2 -1/10], A1)?! not fsolve(@myfun,p0)!

May you please send the code heare that I see, It seems that your code is different to mine.

Thank you

[Walter Roberson]

fsolve(@(p) myfun(p, A1), [2800 -1/2 -1/10])

[azam ghamari]

Walter

it is my whole code :

a myfun mfile:

function F = myfun(p)
F=A1'*P+P*A1+eye(2,2)=0;
p0=[28000 -0.5;-0.5 -0.1];
end

and below code:

now also I write what you write, it gives me error. in past , I wrote "[p,fval]=fsolve(@myfun,p0)" also it gave me the error. By the way I write like this p0=[28000 -0.5;-0.5 -0.1] for the initial condition.May you please put your code here to see it.

Thanks

clear all;
clc;
close all;
A1=[3.4*10^-15 7.6*10^-6;1 4];
%A1=[3.4*10^-15 7.6*10^-6 0 0;-4*10^-16 -6.5*10^-5 0 0;2.3*10^-16 0.00059 -0.1426 0;-8.159*10^-14 -0.05 18.21 0.08];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.977*10^-7 -1.008*10^-5 0 0;9.24*10^-8 0.00045 -0.127 0;-8.55*10^-5 -0.10 65.217 -0.0833];
syms p(1)
syms p(2)
syms p(3)
syms p4
syms p5
syms p6
syms p7
syms p8
syms p9
syms p10
syms P
%P=[p1 p2 p3 p4;p2 p5 p6 p7;p3 p6 p8 p9;p4 p7 p9 p10];
p=[p(1) p(2);p(2) p(3)];
p(1)>0;
det(P)>0;
[p,fval]=fsolve(@(p)myfun(p,A1),p0)

[Walter Roberson]

I already posted the code for myfun in https://www.mathworks.com/matlabcentral/answers/455900-i-want-to-solve-a-10-equation-10-states-with-fsolve-but-it-does-not-give-me-the-answer#answer_370298

[azam ghamari]

I see your code, and when I wanted to exactly write it, it gives me just this:

p1 =

p1

and when I write

function F = myfun(p,A1)

P=[p(1) p(2);p(2) p(3)];

F=A1'*P+P*A1+eye(2,2)=0;

p0=[28000 -0.5;-0.5 -0.1];

end

clear all;
clc;
close all;
A1=[3.4*10^-15 7.6*10^-6;1 4];
%A1=[3.4*10^-15 7.6*10^-6 0 0;-4*10^-16 -6.5*10^-5 0 0;2.3*10^-16 0.00059 -0.1426 0;-8.159*10^-14 -0.05 18.21 0.08];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.977*10^-7 -1.008*10^-5 0 0;9.24*10^-8 0.00045 -0.127 0;-8.55*10^-5 -0.10 65.217 -0.0833];
syms p1 p2 p3 P
sol=fsolve(@(p) myfun(p,A1),p0);
P=vpa([sol.p1,sol.p2,sol.p3]);
display (P)

it also gives me the same

p1=

p1

!! no answer.

however, you write once with solve and once with fsolve. finally which of them?

Thanks

[Walter Roberson]

file myfun.m :

function F = myfun(p, A1)
P=[p(1) p(2);p(2) p(3)];
F=A1'*P+P*A1+eye(2,2);
end

file solve22.m:

A1=[3.4*10^-15 7.6*10^-6;1 4];
p0 = [28000 -1/2 -1/10];
[P,fval]=fsolve(@(P) myfun(P,A1),p0)

file solve44.m:

A1=[3.4*10^-15 7.6*10^-6;1 4];

P = sym('p', [1 3]);

sol = solve(myfun(P,A1));
p = vpa([sol.p1, sol.p2, sol.p3]);
disp(p)

[azam ghamari]

Walter thanks , but none of your codes give me the solution, ther just give me

p1=

p1

. I don't know why!!!!

what is your version of matlab?

B.R

[azam ghamari]

Dear Mr Walter

I wrote your codes(Fsolve), but it gives me this result:

"No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function

values is not near zero, as measured by the default value of the function tolerance.

<stopping criteria details>"

and when I change the initial condition, it does not go and search, it just stops at the same point.

However, fsolve also always gives me some errors that p i snot defined.

[Walter Roberson]

It is going to give you the "ineffective" message if your starting point is not quite close to the actual answer. 28000 might not be close enough: you might have to use 27900 .

When it does work, fsolve() will probably give a result close to the input values instead of giving an exact value.

In other words, fsolve() is not a good way of doing these calculations: if you are not close enough to the answer when you start then it will probably not find any solution, and if you are close enough to the answer when you start then it will probably say "That's close enough" unless you pass in options to force it use tighter tolerances.

The symbolic approach on the other hand finds an exact solution quickly.

[azam ghamari]

Dear Mr Walter

I wrote the solve function but it gives me this error:

Error using solve (line 268)

Specify a variable for which you solve.

Error in sample (line 15)

sol=solve(myfun(P,A1));

(I write the myfun in another mfile and copy it below the main program)

clear all;
clc;
close all;
sym p(1)
sym p(2)
sym p(3)
sym P
A1=[3.4*10^-15 7.6*10^-6;1 4];
p0=[1 -0.5;1 -0.1];
 [P,fval]=fsolve(@(P) myfun(P,A1),p0);
sol=solve(myfun(P,A1));
 p=vpa([sol.p(1),sol.p(2),sol.p(3)]);
 display (P);
-------------------------------------
function F = myfun(p,A1)
P=[p(1) p(2);p(2) p(3)];
F=A1'*P+P*A1+eye(2,2);
end
 
 

[Walter Roberson]

When you did the fsolve(), you overwrote the P variable that I had created by

P = sym('p', [1 3]);

Use a different output variable for fsolve(), or at least postpone the fsolve() until after the solve()

Note: those lines with sym p(1) and so on are wrong and will cause an error. You will also find that sol.p does not have 3 entries.

I do not know why you insist on breaking the working code that I give you.

Answer 3

[azam ghamari]

no I write you code exactly and I comment the line fsolve, wrongly it's written here. However I write your code:

with this myfun code:

function F = myfun(p,A1)

P=[p(1) p(2);p(2) p(3)];

F=A1'*P+P*A1+eye(2,2);

end

but it gives me this error:

Reference to non-existent field 'p'.

Error in sample (line 16)

p=vpa([sol.p(1),sol.p(2),sol.p(3)]);

clear all;
clc;
close all;
A1=[3.4*10^-15 7.6*10^-6;1 4];
P=sym('p',[1 3]);
p0=[1 -0.5;1 -0.1];
 %[P,fval]=fsolve(@(P) myfun(P,A1),p0);
sol=solve(myfun(P,A1));
 p=vpa([sol.p(1),sol.p(2),sol.p(3)]);
 display (P);

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/455900-i-want-to-solve-a-10-equation-10-states-with-fsolve-but-it-does-not-give-me-the-answer#comment_693628 and notice I refer to sol.p1 not sol.p(1)

[azam ghamari]

ok, thanks very much, it works.

so this equations is better to solve with solve not fsolve.... so which equation we can solve from fsolve? (if we know the answer that we start from around that work, it is solved, but we don't know the answer at all. it is the problem)

thanks