How can I solve the non linear algebraic equation system with a parameter

I tried as shown below :

function F=root2d(x)

syms s

F(1)=46.263e+5*x(1)+0.288875*x(3)^2*x(2)*sin(x(4))-0.187375*x(2)^2*x(3)*sin(x(6))-1.5e+3*sin(x(5)); F(2)=237.807241*s*x(1)-0.7765*x(1)^3-2.07225*x(1)*x(2)^2-0.288875*x(1)^2*x(2)*cos(x(4))+0.187375*x(2)^2*x(3)*cos(x(6))-4.341625*x(1)*x(3)^2+1.5e3*cos(x(5)); F(3)=249.85e5*x(2)-0.09625*x(1)^3*sin(x(4))+0.374875*x(1)*x(2)*x(3)*sin(x(6)); F(4)=2311.76211*x(2)*(z-801)-10.7815*x(2)^3-2.07225*x(1)^2*x(2)-0.09625*x(1)^3*cos(x(4))+0.374875*x(1)*x(2)*x(3)*cos(x(6))-16.2345*x(3)^2*x(2); F(5)=938.16e5*x(3)+0.187375*x(1)*x(2)^2*sin(x(6)); F(6)=(11252.3953*s+3866749.2347199973)*x(3)+0.187375*x(1)*x(2)^2*cos(x(6))-4.341625*x(1)^2*x(3)-16.2345*x(2)^2*x(3)-51.815125*x(3)^3;

options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt); fun = @root2d; x0=[1,1,1,1,1,1]; x=fsolve(fun,x0,options)

This is the Error I am geting :

The following error occurred converting from sym to double: DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use VPA.

Error in root2d (line 6) F(2)=237.807241*s*x(1)-0.7765*x(1)^3-2.07225*x(1)*x(2)^2-0.288875*x(1)^2*x(2)*cos(x(4))+0.187375*x(2)^2*x(3)*cos(x(6))-4.341625*x(1)*x(3)^2+1.5e3*cos(x(5));

Error in fsolve (line 218) fuser = feval(funfcn{3},x,varargin{:});

Caused by: Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Answer 1

[Star Strider]

You can use the matlabFunction function to create an anonymous function from ‘F’ so you can use it with fsolve. However, you still need to provide values for ‘s’ and ‘z’ in order to solve it for ‘x’ (that appears in the ‘root2d’ anonymous function as ‘in1’).

The Code —

syms s x z 
x = sym('x', [1,6]);
F(1)=46.263e+5*x(1)+0.288875*x(3)^2*x(2)*sin(x(4))-0.187375*x(2)^2*x(3)*sin(x(6))-1.5e+3*sin(x(5)); 
F(2)=237.807241*s*x(1)-0.7765*x(1)^3-2.07225*x(1)*x(2)^2-0.288875*x(1)^2*x(2)*cos(x(4))+0.187375*x(2)^2*x(3)*cos(x(6))-4.341625*x(1)*x(3)^2+1.5e3*cos(x(5)); 
F(3)=249.85e5*x(2)-0.09625*x(1)^3*sin(x(4))+0.374875*x(1)*x(2)*x(3)*sin(x(6)); 
F(4)=2311.76211*x(2)*(z-801)-10.7815*x(2)^3-2.07225*x(1)^2*x(2)-0.09625*x(1)^3*cos(x(4))+0.374875*x(1)*x(2)*x(3)*cos(x(6))-16.2345*x(3)^2*x(2); 
F(5)=938.16e5*x(3)+0.187375*x(1)*x(2)^2*sin(x(6)); 
F(6)=(11252.3953*s+3866749.2347199973)*x(3)+0.187375*x(1)*x(2)^2*cos(x(6))-4.341625*x(1)^2*x(3)-16.2345*x(2)^2*x(3)-51.815125*x(3)^3;
root2d = matlabFunction(F, 'Vars',{[x],s,z})
s = rand;                                                                               % Provide Correct Value (Scalar)
z = rand;                                                                               % Provide Correct Value (Scalar)
options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt); 
fun = @(in1)root2d(in1,s,z); 
x0=[1,1,1,1,1,1]; 
x=fsolve(fun,x0,options)

[Thimé Cianyi]

Thanks dear Strider, I’ve tried to run the above code and it gives only for given value of s. I made a mistake when writing the function F(4), it is s instead of z.

Right now I am concerned in getting the unknown for s in range of -200 and 600. And plotting x(i) vs s, with i=1,...,6.

What could be the code for this ?

[Star Strider]

Try this:

syms s x 
x = sym('x', [1,6]);
F(1)=46.263e+5*x(1)+0.288875*x(3)^2*x(2)*sin(x(4))-0.187375*x(2)^2*x(3)*sin(x(6))-1.5e+3*sin(x(5)); 
F(2)=237.807241*s*x(1)-0.7765*x(1)^3-2.07225*x(1)*x(2)^2-0.288875*x(1)^2*x(2)*cos(x(4))+0.187375*x(2)^2*x(3)*cos(x(6))-4.341625*x(1)*x(3)^2+1.5e3*cos(x(5)); 
F(3)=249.85e5*x(2)-0.09625*x(1)^3*sin(x(4))+0.374875*x(1)*x(2)*x(3)*sin(x(6)); 
F(4)=2311.76211*x(2)*(s-801)-10.7815*x(2)^3-2.07225*x(1)^2*x(2)-0.09625*x(1)^3*cos(x(4))+0.374875*x(1)*x(2)*x(3)*cos(x(6))-16.2345*x(3)^2*x(2); 
F(5)=938.16e5*x(3)+0.187375*x(1)*x(2)^2*sin(x(6)); 
F(6)=(11252.3953*s+3866749.2347199973)*x(3)+0.187375*x(1)*x(2)^2*cos(x(6))-4.341625*x(1)^2*x(3)-16.2345*x(2)^2*x(3)-51.815125*x(3)^3;
root2d = matlabFunction(F, 'Vars',{[x],s})
sv = linspace(-200, 600, 9);                                                                % Create ‘s’ Vector
options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt); 
for k = 1:numel(sv)
fun = @(in1)root2d(in1,sv(k)); 
x0=[1,1,1,1,1,1]; 
xm{k} = fsolve(fun,x0,options);
end

figure
for k = 1:numel(sv)
subplot(numel(sv),1,k)
plot((1:6), xm{k})
title(sprintf('s = %d', sv(k)))
set(gca, 'XTickLabel', (1:6))
end
fpos = get(gcf, 'Position');
set(gcf, 'Position', fpos+[0 -500 0 +500])                                                  % Set Figure Window To Show Subplots Easily

[Thimé Cianyi]

Thanks dear Strider, From this stage the output looks like below :

But what I really want is plotting x(i, i=1,2,3,4,5,6) versus ‘s’ all the x(i) within the same figure but different colors.

From where I should modify the code?

Thank you.

[Star Strider]

The loop needs to change from saving a cell array to a numeric array, so it becomes:

for k = 1:numel(sv)
    fun = @(in1)root2d(in1,sv(k));
    x0=[1,1,1,1,1,1];
    xm(k,:) = fsolve(fun,x0,options);
end

and the plot call becomes:

figure
plot((1:6), xm)
grid
lgdc = sprintfc('s = %d', sv);
legend(lgdc, 'Location','SW')

I used the subplot approach because the magnitudes of the ‘x’ values change significantly, and it shows the full range of values for all of them. I initially saved ‘xm’ as a cell array because I did not know if the outputs would change. They do not, so a numeric array is appropriate.

[Thimé Cianyi]

Thank you, but still no plotting in ordinates x_i and in abscisse s

I’ve tried to modify the code like this

root2d= matlabFunction(F,'vars',{[x],s}); sv = linspace(-0.02,0.06); options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt); for k = 1:numel(sv) fun = @(in1)root2d(in1,sv(k)); x0 = [1,1,1,1,1,1]; xm(k,:) =fsolve(fun,x0,options); end

figure plot(sv,xm) grid lgdc = sprintfc('x= %d', xm); legend(lgdc, 'location', 'SW')

And I’ve got 6 curves each one represents x_i for -0.02<= s <=0.06

But I didn’t manage to label all the curve as x_i, how can I do this ?

Thank you again

[Star Strider]

I am not certain what you want, since you did not mention how you wanted the values to be plotted.

Note that:

sv = linspace(-0.02,0.06);

is going to create 100 values for ‘sv’, so that is going to be extremely difficult to fit in a legend entry.

Try this:

sv = linspace(-0.02, 0.06);
for k = 1:numel(sv)
    fun = @(in1)root2d(in1,sv(k));
    x0=[1,1,1,1,1,1];
    xm(k,:) = fsolve(fun,x0);
end

figure
plot(xm, sv)
grid
xlabel('x')
ylabel('s')

I removed the ‘options’ that produced the fsolve plots because it was taking too long with 100 iterations.

[Thimé Cianyi]

Thanks I think this is working and the output is really what I want when changing plot(sv, am).

An extra thing I want to know is : Suppose I have three values for F1 : 30, 50 and 90. And I want to plot a_1 versus s knowing that s goes from -0.02 to 0.06. I want the output to be three curves of a_1, each curve with his value of F1. What could be the last part of the code ?

[Star Strider]

i do not understand what you want to do.

What is the relation between ‘a_1’ (that just now appeared for the first time), ‘F(1)’, and ‘s’?

[Thimé Cianyi]

I am sorry I thought my code was uploaded here when asking question.

This is the code :

syms sigma_2 F_1 omega_1 = 15.421; omega_2 = 49.970; omega_3 = 104.240; sigma_12 = 801; sigma_13 = - 964; c_1 = 100; c_2 = 100; c_3 = 100; F_1 = 30; q_1 = -6.212; q_2 = -16.578; q_3 = -2.311; q_4 = -86.252; q_5 = -16.578; q_6 = -0.770; q_7 =1.499; q_8 = 2.999; q_9 = 1.499; q_10 = -34.733; q_11 = -129.876; q_12 = -34.733; q_13 = -129.876; q_14 = -414.521; a= sym('a',[1,6]);

F(1) = omega_1 * c_1 * a(1) - (q_3 * (a(3))^2 * a(2) * sin(a(4)) + q_7 * (a(2))^2 * a(3) * sin(a(6)))/8 - (force * sin(a(5)))/2; F(2) = (omega_1)^2 * sigma_2 * a(1) + (q_1 * (a(1))^3 + q_2 * a(1) * (a(2))^2 + q_3 * (a(1))^2 * a(2) * cos(a(4)) + q_7 * (a(2))^2 * a(3) * cos(a(6)) + q_10 * a(1) * (a(3))^2)/8 + (force * cos(a(5)))/2; F(3) = omega_2 * c_2 * a(2) + (q_6 * (a(1))^3 * sin(a(4)) - q_8 * a(1) * a(2) * a(3) * sin(a(6)))/8; F(4) = 3* omega_1 * omega_2 * a(2) * (sigma_2 - sigma_12) + (q_4 * (a(2))^3 + q_5 * (a(1))^2 * a(2) + q_6 * (a(1))^3 * cos(a(4)) + q_8 * a(1) * a(2) * a(3) * cos(a(6)) + q_11 * (a(3))^2 * a(2))/8; F(5) = omega_3 * c_3 * a(3) + (q_9 * a(1) * (a(2))^2 * sin(a(6)))/8 ; F(6) = omega_1 * omega_3 * a(3) * (7*sigma_2 - 6*sigma_12 - ((omega_1 + 2 * omega_2)* sigma_13 / omega_1)) + (q_9 * a(1) * (a(2))^2 * cos(a(6)) + q_12 * (a(1))^2 * a(3) + q_13 * a(3) * (a(2))^2 + q_14 * (a(3))^3)/8;

root2d= matlabFunction(F,'vars',{[a],sigma_2}); sigma2v = linspace(-200,600); %forcev = [3e+5 5e+5 9e+5]; options = optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt);

for k = 1:numel(sigma2v) fun = @(in1)root2d(in1,sigma2v(k)); a0 = [1,1,1,1,1,1]; am(k,:) =abs(fsolve(fun,a0,options)); end figure

plot(sigma2v,am(:,1)) grid

xlabel('sigma_2'); ylabel('Amplitude a'); ———————-

F1 and s are parameters in the equation system, before I was plotting a versus s with F1=30 only, but now I want to take into account that F1 vary from 30, 50 and 90. What I want is to plot one unknown a_i versus s with différents values of F1 in order to get three curves on the same figure.

[Star Strider]

I cannot read the code you posted. You need to format it as code.

This also does not appear to be anything similar to the code you originally posted.

‘F1 and s are parameters in the equation system ...’

Not the one you originally posted. The only ‘F’ I see are the elements of the vector of equations you are solving with fsolve.

[Thimé Cianyi]

Okay I am taking a clear example. Let’s consider the system of equations with 2 parameters s and m :

F(1) = 2s x(1) + 4m x(2)^2 F(2) = sin(x(1)) + exp(1.8m x(2))

And I want to solve it using fsolver and plotting the solution x(1) and (2) versus s with -2<s<5 x(1) for each value of m ( m=1,2,3).

The code should be able to plot x(1) versus s with -2<s<5 and for each value of m ( m=1,2,3). And give me the possibility of modifying it to plot the same for x(2).

From this example I can be able to solve the complicated one that have.

I hope this to be more clear, thank you.

[Star Strider]

Try this:

F = @(x,s,m) [2*s.*x(1) + 4*m.*x(2).^2; sin(x(1)) + exp(1.8*m.*x(2))];
mv = 1:3;
sv = -2:5;
for k1 = 1:numel(mv)
    for k2 = 1:numel(sv)
        xm(k1,k2,:) = fsolve(@(x)F(x,sv(k2),mv(k1)), [1 1]);
    end
end

figure
plot(sv, xm(:,:,1), '--')
hold on
plot(sv, xm(:,:,2), '.-')
hold off
xlabel('s')
ylabel('x')
x1sc = sprintfc('x_1, m = %d',mv);
x2sc = sprintfc('x_2, m = %d',mv);
legend([x1sc, x2sc], 'Location','SE')

[Thimé Cianyi]

Thank you dear Star, I’ve got it and I implemented the idea on my complicated code one, and this is in if the output:!

[Star Strider]

As always, my pleasure!

I am quite happy that we were able to solve this!

Answer 2

[John D'Errico]

You CANNOT use fsolve with a symbolic parameter in the problem. CANNOT. Fsolve is a numerical tool. It finds a root (if it can that is purely numerical.)

If you are willing to choose some specific value of s, then the result will be 6 equations in 6 unknowns. There will usually be multiple solutions to such a problem, especially when you have trig functions in there. Fsolve might find ONE of them, IF it converges to a valid solution. The result you get will be dependent on the starting values you provide.

You have a system of 6 equations in 6 unknowns, with one parameter that you want to leave in there, so the solutinos will be a function of s. Unfortunately, your system is complicated enough that it will surely not have any analytical solutino at all. Powers of the x(i), as well as having those variables inside and outside of trig functions mean that there will never be an analytical solution. So even if you wanted to try to use solve, where the result could be a function of a parameter, it will fail.