Plotting the results of a nonlinear system of equations as a parameter changes

Hey guys!

I,m quite new to matlab.

I have a system of four nonlinear equations for which I want to plot the results as a function of a variable. Can anybody help me on that?

function F = myfun(x)
w1 = 1;
w2 = 3;
a1 = -24.851;
a3 = -66.319;
a2 = -9.247;
a8 = -345.019;
a4 = 0;
a6 = -66.319;
a5 = -3.082;
a7 = 0;
mu1 = 0;
mu2 = 0;
f1 = 0;
f2 = 0;
S1 = 1;
S2 = 0.01;
F = [ 8 * w1 * mu1 * x(1) - a2 * x(1)^2 * x(2) * sin(x(3)) - 4 * f1 * sin(x(4));
    8 * w2 * mu2 * x(2) + a5 * x(1)^3 * sin(x(3));
    8 * w1 * x(1) * S2 + (3* a1 * x(1)^2 + 2 * a3 * x(2)^2) * x(1) + a2 * x(1)^2 * x(2) * x(2) * cos(x(3)) + 4* f1 * cos(x(4));
    8 * w2 * x(2) * (3 * S2 - S1) + (3 * a8 * x(2)^2 + 2 * a6 * x(1)^2) * x(2) + a5 * x(1)^3 * cos(x(3))];
end

Here is my code. the unknowns are x(1),x(2),x(3) and x(4). What I want is:

1. Solve the nonlinear system and show the answer.
2. Plot x(1) vs S1 parameter as S1 changes from -50 to 50.

Is this possible to do in matlab?

Thanks in advance.

Answer 1

[John D'Errico]

Well, you did get to a start, and this does not look directly like homework.

I'll rewrite your function, making it a function of only S1, then call solve inside.

function res = myfun(S1)
syms x1 x2 x3 x4
w1 = 1;
w2 = 3;
a1 = -24.851;
a3 = -66.319;
a2 = -9.247;
a8 = -345.019;
a4 = 0;
a6 = -66.319;
a5 = -3.082;
a7 = 0;
mu1 = 0;
mu2 = 0;
f1 = 0;
f2 = 0;
S2 = 0.01;
E1 = 8 * w1 * mu1 * x1 - a2 * x1^2 * x2 * sin(x3) - 4 * f1 * sin(x4);
E2 = 8 * w2 * mu2 * x2 + a5 * x1^3 * sin(x3);
E3 = 8 * w1 * x1 * S2 + (3* a1 * x1^2 + 2 * a3 * x2^2) * x1 + a2 * x1^2 * x2 * x2 * cos(x3) + 4* f1 * cos(x4);
E4 = 8 * w2 * x2 * (3 * S2 - S1) + (3 * a8 * x2^2 + 2 * a6 * x1^2) * x2 + a5 * x1^3 * cos(x3);
res = solve(E1,E2,E3,E4,x1,x2,x3,x4);

It looks like there are generally 9 independent solutions, 3 of them are not a function of S1 at all. As well, x(2) always seems to come out as zero in all solutions.

Lets take a look at a few cases for different values of S1.

res = myfun(-50);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[      0, 0,         0,    -1.0771]
[      0, 0,         0,     1.0771]
[      0, 0, -0.032758, -9.0236e-8]
[      0, 0,         0,          0]
[ 3.1416, 0,         0,    -1.0771]
[ 3.1416, 0,         0,     1.0771]
[ 3.1416, 0,         0,          0]
[      0, 0,    1.5708,          0]
[ 3.1416, 0,  0.032758, -9.0236e-8]

res = myfun(-30);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[      0, 0,         0,   -0.83445]
[      0, 0,         0,    0.83445]
[      0, 0,         0,          0]
[ 3.1416, 0,  0.032758, -1.5035e-7]
[ 3.1416, 0,         0,   -0.83445]
[ 3.1416, 0,         0,    0.83445]
[      0, 0, -0.032758, -1.5035e-7]
[ 3.1416, 0,         0,          0]
[      0, 0,    1.5708,          0]

res = myfun(0);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[      0, 0,         0,   -0.026374]
[      0, 0,         0,    0.026374]
[ 3.1416, 0,  0.032757, -0.00018753]
[      0, 0,         0,           0]
[      0, 0, -0.032757, -0.00018753]
[ 3.1416, 0,         0,   -0.026374]
[ 3.1416, 0,         0,    0.026374]
[ 3.1416, 0,         0,           0]
[      0, 0,    1.5708,           0]

res = myfun(25);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[      0, 0,  0.032758, -1.8073e-7]
[      0, 0,         0,  -0.76091i]
[      0, 0,         0,   0.76091i]
[      0, 0,         0,          0]
[ 3.1416, 0,         0,          0]
[      0, 0,    1.5708,          0]
[ 3.1416, 0, -0.032758, -1.8073e-7]
[ 3.1416, 0,         0,  -0.76091i]
[ 3.1416, 0,         0,   0.76091i]

res = myfun(10);
vpa([res.x1,res.x2,res.x3,res.x4],5)
ans =
[      0, 0,         0,  -0.48081i]
[      0, 0,         0,   0.48081i]
[      0, 0,         0,          0]
[ 3.1416, 0,         0,          0]
[      0, 0,    1.5708,          0]
[ 3.1416, 0, -0.032758, -4.5249e-7]
[      0, 0,  0.032758, -4.5249e-7]
[ 3.1416, 0,         0,  -0.48081i]
[ 3.1416, 0,         0,   0.48081i]

Be careful, because many of the solutions are the same, just reported in a different sequence.

So x2 is always zero for all the cases I've tried. You could probably show that to be necessary by looking at the equations.

Those 3.1416 values are pi, rounded off to 5 digits by vpa. In fact, the symbolic toolbox shows it to be exactly pi. And pi/2 is 1.5708.

So some solutions in common to all are [0,0,0,0], [0,0,pi/2,0], and [pi,0,0,0]. I'd want to ignore all of those apparently trivial solutions.

Next, while all solutions seem to be real for negative S1, at some point it seems that solutions start turning imaginary for at least some values of positive S1.

There appear to be 4 solutions of the form

[0,0,0, x(4)]
[0,0,0,-x(4)]
[pi,0,0, x(4)]
[pi,0,0,-x(4)]

There x(4) is a function of S1. With a little algebra, that relationship could surely be uncovered.

Finally, there always seem to be a pair of solutions of the form

[ 0,0,-x3,x4]
[pi,0, x3,x4]

Here x3 and x4 are both functions of S1. Again, some simple algebra should recover that relationship. That seems to account for all 9 basic solutions.

Quick inspection of your coefficients has many of them as zero, f1, for example. A result of that is Equation 2:

E2
E2 =
-(1541*x1^3*sin(x3))/500

E2 is zero ONLY when x1 is zero, OR sin(x3) is zero. That happens at x1==0, x3==0, or x3==k*pi, for any integer k.

You can do the rest.

[SorSep]

Thank you very much sir. Now is there a way to plot x(1) versus S1 for this set of equations?

[John D'Errico]

Reread what I wrote. X(1) is NEVER a function of S1. NEVER. In all solutions, there is no relationship between S1 and x(1).

So, I suppose if you really, desperately, absolutely want to plot a horizontal line, feel free to do so. :)

You can spend a few moments doing some more analysis to derive the relationship between S1 and X3 & X4.

[SorSep]

Thanks again sir. but what you're proposing feels a kinda incorrect...

These are the four equations. The unknowns are a1,a2,gamma1 and gamma2. I want to plot a1 and a2 as a function of sigma2. The final plot should be something like this:

as you see, a1 IS actually a function of sigma2. All is want is to be able to plot a1 as a function of sigma2 as below. The values for the parameters are the same as what I mentioned above. sigma1 might have a different value than 1 though.

[SorSep]

I believe I made a mistake sir. S2 should not have a numeric value abviously. So my four equations would be like:

function F = myfun(x)
w1 = 1;
w2 = 3;
a1 = -24.851;
a3 = -66.319;
a2 = -9.247;
a8 = -345.019;
a4 = 0;
a6 = -66.319;
a5 = -3.082;
a7 = 0;
mu1 = 0;
mu2 = 0;
f1 = 0;
f2 = 0;
S1 = 1;
F = [ 8 * w1 * mu1 * x(1) - a2 * x(1)^2 * x(2) * sin(x(3)) - 4 * f1 * sin(x(4));
    8 * w2 * mu2 * x(2) + a5 * x(1)^3 * sin(x(3));
    8 * w1 * x(1) * S2 + (3* a1 * x(1)^2 + 2 * a3 * x(2)^2) * x(1) + a2 * x(1)^2 * x(2) * x(2) * cos(x(3)) + 4* f1 * cos(x(4));
    8 * w2 * x(2) * (3 * S2 - S1) + (3 * a8 * x(2)^2 + 2 * a6 * x(1)^2) * x(2) + a5 * x(1)^3 * cos(x(3))];
end

[SorSep]

and just to clarify. a1 , a2 ,...,a8 in the code actually correspond to alpha1, alpha2 ,...,alpha8 and w1 and w2 are omega1 and omega2 respectively.