Fsolve with Loop and to store variable

I am having the following set of 4 non linrear equations with 6 variable x1,x2,x3,x4,x5,x6

x1-(5*10^3)*x2-353=0
2*x2*(15*10^-3)-x3-x4-x5=x1=0
4000*x3*x2-x6^2=0
8.2*x6+x4*x2-1259=0

Initial guess is x1=603,x2=5*10^4,x3=9.5*10^-6,x4=0,x5=1374,x6=1373

For the equations to solve I have used the following code :

x0 = [603 5*10^4 9.5*10^-6 0 1374 1373]
F = @(x) [x(1)-(5*10^3*x(2))-353;
 0.03*x(2)-(2*x(2)*x(3))-(2*x(2)*x(4))-(2*x(2)*x(5))+(2*x(2)*(x1));
4000*x(3)*x(2)-x(6)^2;
8.2*x(6)+x(4)*x(2)-1259];
x0 = [603 5*(10^4) 9.5*(10^-6) 0 1374 1373];
options = optimoptions('fsolve','Display','iter');
[x,fval] = fsolve(F,x0,options)

error is unrecognised x1...............

Further I want a loop to be incorporated:

This equations has to ve solved iteratively for (x4)i+1= (x4)i+2.2*10^-6*(delta T) where T = [0.1:1:20]

I am interested to store and plot x1,x2,x3,x4,x5,x6 vs T

Any help will be very much appreciated.... I tried to intiate the solution with fsolve ... but I am not in position to call functions, execute the for loop and store the iteration results for plotting.

Thanking you in advance !

[Walter Roberson]

2*x2*(15*10^-3)-x3-x4-x5=x1=0

please recheck that. It implies that x1 must always be 0.

[Walter Roberson]

 0.03*x(2)-(2*x(2)*x(3))-(2*x(2)*x(4))-(2*x(2)*x(5))+(2*x(2)*(x1));

is a fair bit different than the second equation you posted.

[Walter Roberson]

And your x1 should be x(1)

Answer 1

[Walter Roberson]

You are required to iterate through given x(4) values. So start solving the equations in terms of x(4)

syms x1 x2 x3 x4 x5 x6
eqn = [
    x1-(5*10^3)*x2-353==0
    2*x2*(15*10^-3)-x3-x4-x5+x1==0   %equation was corrupt, had to guess what it was
    4000*x3*x2-x6^2==0
    8.2*x6+x4*x2-1259==0]

eqn = 

sol = solve(eqn, [x1, x2, x3, x5]);
SOL = simplify([sol.x1; sol.x2; sol.x3; sol.x5])

SOL = 

residue = sum((SOL - [603; 5*10^4; 9.5*10^-6; 1374]).^2)

residue = 

%do not display, it is long long equations
bestx6 = solve(diff(residue, x6),x6, 'maxdegree', 4);
x4vals = linspace(1e-8, 2.2E-6*20, 100);
X6 = double(subs(bestx6(1:2), x4, x4vals)).';  %should be 100 x 2
subplot(2,1,1); plot(x4vals, real(X6(:,1:2))); title('real'); legend({'1', '2'})
x4vals = logspace(-20, -8, 50);
X6 = double(subs(bestx6(1:2), x4, x4vals)).';  %should be 100 x 2
subplot(2,1,2); plot(x4vals, real(X6(:,1:2))); title('real'); legend({'1', '2'})

So what is the best solution, the one in which the all of the variables are as close as possible to the starting values? Answer: it is the second out of four x6 solutions, evaluated at the largest x4 value. And to the accuracy you would normally use for double precision, you cannot tell the solutions apart, so it would be valid to use either one.

What did we calculate here? We found formulas for x1, x2, x3, x5 in terms of x4 and x6, x4 is given, which leaves x6. So which x6 value gives us the closest possible result to the initial values given as part of the problem? We can optimize for least squared error using calculus techniques that give us optimal x6 in terms of x4.

Then to find the overall best, we want the x6 closest to the initial value given for x6. But it turns out that the optimal x6 values are a long way from the initial value given for x6 and are to be found over a very narrow range considering the range of x4 values that are to be examined.

Since you have 4 equations in 6 unknowns, one of which is fixed at any one point, we need to give meaning to assigning the "best" value to the remaining variable.

There are other ways that you could measure "best" solution -- for example you might find some way to fit the value of x6 relative to the initial value (1593) in to the measure.

[Jhon Paul]

Thank you very much for your prompt reply.. I guess I have made some mistakes in posting the equations... will follow the steps given by you and update.

Lastly .. I have been following all your response .. A true Gem ! :)

[Jhon Paul]

Dear Walter,

I went through the code and it really works fine .... However after getting the solutions of each variable, you are directing it towards optimisation. As I see from the solution, all the solution values are close and my intention is not to optmise the problem. So I can take any one solution as you have suggested. I have corrected the code, which you rightly have pointed out. Is there any way out that now I incorporate a loop with one of the solution of x4. My intention is to plot all the variable Xi , while the value of x4 is incremented in steps of (x4)i+1= (x4)i+ 5*10^-3.

syms x1 x2 x3 x4 x5 x6
Eqn = [
    x1-5000*x2-353==0
    x1+(3/100)*x2-2*x2*x3-2*x2*x4-x5==0   
    4000*x2*x3-x6^2==0
    8.2*x6+x4*x2-11259==0]
sol = solve(Eqn, [x1, x2, x3, x6]);
SOL = simplify([sol.x1; sol.x2; sol.x3; sol.x6])

[Walter Roberson]

syms x1 x2 x3 x4 x5 x6
Eqn = [
    x1-5000*x2-353==0
    x1+(3/100)*x2-2*x2*x3-2*x2*x4-x5==0   
    4000*x2*x3-x6^2==0
    8.2*x6+x4*x2-11259==0]

Eqn = 

sol = solve(Eqn, [x1, x2, x3, x6]);
SOL = simplify([sol.x1, sol.x2, sol.x3, sol.x6])

SOL = 

x5_bound = solve(56157500*x4^2 - 500*x4^2*x5-308054348315*x4+420255043015129, x5)

x5_bound = 

x4_bound = solve(x5_bound)

x4_bound = 

vpa(x4_bound)

ans = 

S5 = subs([sol.x1, sol.x2, sol.x3, [x4;x4], [x5;x5], sol.x6], x5, x5_bound);
t1 = vpa(subs(S5,x4,2500))

t1 = 

t2 = vpa(subs(S5, x4, 2750))

t2 = 

t3 = vpa(subs(S5, x4, 3000))

t3 = 

[Walter Roberson]

What the above tells you is that when you solve for the variables in terms of x4 and x5, then there is an upper bound for x5 in terms of x4, x5_bound, with large x5 values leading to complex results.

The t1, t2, t3 show you that there are a wide range of possible results for x1, x2, x3, x6, depending on the x4 and x5 values.

[Walter Roberson]

my intention is not to optmise the problem.

That is because you do not yet understand what is being done.

The original question asks you to solve the equations iteratively given initial values. You have four equations in six variables, one of which (x4) will be constant at any one time. So at any one time, you have four equations in five free variables. In such a situation, you can solve for any four of the variables in terms of the firth variable.

Now, because you are given initial values, you need to figure out which value of the fifth variable will give you a result that is closest to the initial values that you are given. And that is an optimization problem. I showed you steps you can use (with calculus) to optimize the remaining variable to get solutions as close as possible to the initial values.

The initial values given are

[603, 5E4, 9.5E-6, 0, 1374, 1373]

In my posting two above, I show that with a particular choice of x5 in terms of x4 (the largest x5 consistent with x4), that depending on the x4, you can easily get x1 results anywhere from -18236 to 22872. Considering the initial value of x1 is 603, are those all of equal value to you?

[Walter Roberson]

syms x1 x2 x3 x4 x5 x6
eqn =[
    x1-5000*x2-353==0
    x1+(3/100)*x2-2*x2*x3-2*x2*x4-x5==0   
    4000*x2*x3-x6^2==0
    8.2*x6+x4*x2-11259==0]

eqn = 

sol = solve(eqn, [x1, x2, x3, x6]);
SOL = simplify([sol.x1, sol.x2, sol.x3, [x4;x4], [x5;x5], sol.x6])

SOL = 

x0 = [603, 5*(10^4), 9.5*(10^-6), 0, 1374, 1373];
residue1 = expand(sum((SOL(1,:) - x0).^2))

residue1 = 

[n1, d1] = numden(residue1);
vpa(n1)

ans = 

residue2 = expand(sum((SOL(2,:) - x0).^2))

residue2 = 

[n2, d2] = numden(residue2);
vpa(n2)

ans = 

%do not display, it is long long equations
bestx51 = solve(diff(n1, x5), x5);
bestx52 = solve(diff(n2, x5), x5);
x4vals = linspace(5e-6, 2.2E-6*20, 100);
X51 = double(subs(bestx51(1:2), x4, x4vals)).';  %should be 100 x 2
subplot(4,1,1); plot(x4vals, X51(:,1)); title('branch 1A');
subplot(4,1,2); plot(x4vals, X51(:,2)); title('branch 1B');
X52 = double(subs(bestx52(1:2), x4, x4vals)).';  %should be 100 x 2
subplot(4,1,3); plot(x4vals, X52(:,1)); title('branch 2A');
subplot(4,1,4); plot(x4vals, X52(:,2)); title('branch 2B')

When solving for the various x variables, there are two branches -- essentially you are solving a quadratic in terms of x4 at that point.

Then later you get to the point where you want to find the optimal x5 value to keep the distance from the initial points as small as possible. You are solving x5 in terms of x4, and it is a quartic (degree 4). But two of the roots are imaginary valued, so you have two real values for x5 in terms of x4; and that is still within the context that there are two solutions of the equations in terms of x4.

Plotting (the branch plots) shows that for each of the two overall branches, there is one branch where the optimal x5 (in terms of x4) is nearly constant, approximately -591. That is not especially close to the 1373 initial value originally given, but neither is it a huge difference.

Plotting also shows that for each of the two overall branches, there is one branch where the optimal x5 (in terms of x4) is very dependent on x4; indeed for x4 near 0, x5 increases arbitrarily, as there are divisions by x4^2

It is possible, and would need further research, that for each of the two major solutions sets, that one of the two real branches corresponds to the minima (best x5) and that the other corresponds to a maxima (worst x5). Which is which might differ between the two major solution sets.

[Jhon Paul]

Thank you so much .. I will go through the details carefully ... and update. Thanks for bearing my non-sense.

[Walter Roberson]

It is likely that the expectation is that the system would be solved numerically. However because there are four equations and five free variables, the numeric solutions will depend heavily upon what initial conditions you use and which algorithm you use.

At that point you would need some way of deciding which set of numeric solutions was the "best" solution out of the infinite number of possible solutions. Which, one way or another, becomes an optimization problem.

[Jhon Paul]

Dear Walter Sir,

I totally agree with you and incidentally , the problem has an advantage that x6 almost remain close to the initial value sol.x6 . This means as you have tried to explain , that all the variables can be written in terms of x4 and x6.

I am trying to understand the solution provided by you ... an easy way to incportate the optimisation part through calculus... by minimsing the residuals. I think its almost done in the solution you have given.

Only now, the solution of each variable to eb evaluated with a loop :

Initialisation : sol.x4 = 0

time = 1 to 20E5 (divided into 10 intervals)

velocity= 0.3

Increment : x4vals= sol.x4+ 0.3*linspace(1,20E5,10)

Plot (syms xi,t) where i = 1 to 6

That will do the whole job !

Thanks in advance.

[Walter Roberson]

20E5?? The original request involved 2.2E-6 * 20 which is a maximum of 4.4E-5 and your new range is 5E10 times larger.

[Jhon Paul]

Sorry a correction in the increment part -- you are right !

Increment : x4vals= sol.x4+ (2.2X10^-6)*linspace(1,20E5,10) %% 0.3 was wrong%%

and I can run the whole thing either for time 20E5 or 20 secs. So basically the physics of the problem should not have any change with time . May be you are concerned of the computation time.

I explain you this ...2.2X10-6 refers to the velocity in m/sec. So basically the velocity starts at 0.02 mm/sec. Then when we multiply by time which can be either 20 secs ( as given earlier ) or 900 secs which is the time scale of the problem.

In other words :

Increment : x4vals= sol.x4+ (2.2X10^-6)*linspace(1,900,10)

[Walter Roberson]

Okay, so where is linspace(1,20E5,10) coming from? If 900 seconds is linspace(1,900,10) then 20 seconds would be linspace(1,20,10)

[Jhon Paul]

As I said I can solve the problem either in 20 secs or 900 secs . That won’t change the solution pattern, as the physics remains the same .

[Jhon Paul]

Shall I also try to optimise x4 using optimisation toolbox. I will also try to make the problem dimensionless. I think that will remove all ambiguity and will help to run the problem for a short timescale .

[Walter Roberson]

No, do not optimize x4. The problem statement requires that you use

T = 0.1:1:20;
x0 = [603, 5*(10^4), 9.5*(10^-6), 0, 1374, 1373];
x4_values = x0(4) + 2.2*10^-6 * T;
for K = 1 : length(x4_values)
    x4 = x4_values(K);
    compute with this x4
end

[Jhon Paul]

When I was trying to see the sentiveness of the solution X51 and X52, I observed that both of them are insensitive to the input parameters. Also if I change the x4 range to x4vals = linspace(0.1, 7.83E-3*20, 100);then also I obtain the same branch1A,2A solution.

Does it shows that something I missed as input?

[Jhon Paul]

syms x1 x2 x3 x4 x5 x6
eqn = [
    200*x1-x2-70600 ==0
    x1+3*10^-2*x2*x3-2*x2*x4-x5==0   
    x2*x3-2.5*10^-4*x6^2==0
    x6+0.12*x2*x4-1373==0]
sol = solve(eqn, [x1, x2, x3, x5]);
SOL = simplify([sol.x1; sol.x2; sol.x3; sol.x5])
residue = sum((SOL - [603; 5*10^4; 9.5*10^-3; 1153]).^2)
bestx6 = solve(diff(residue, x6),x6, 'maxdegree', 4)
x4vals = linspace(1e-8, 15E-3, 100)
X6 = double(subs(bestx6(1:6), x4, x4vals)).';  %should be 100 x 2
subplot(3,1,1); plot(x4vals, real(X6(:,1:2))); title('real'); legend({'1', '2',})
subplot(3,1,2); plot(x4vals, real(X6(:,3:4))); title('real'); legend({'3', '4',})
subplot(3,1,3); plot(x4vals, real(X6(:,5:6))); title('real'); legend({'5', '6',})
%% 
T = 0.1:1000:18000;
x0 = [603, 5*(10^4), 9.5*(10^-3), 0, 1153, 1373];
x4_values = x0(4) + 2.2*10^-6 * T;
for K = 1 : length(x4_values)
    x4 = x4_values(K);
end 
plot (T,x4_values)

Dear Walter....

I got it solved .... as you suggested .. there is this one solution of 6 which is bestx6(:,2). Can you now please help me to get the reamining part... x1 to x6..

Also please help to get the expression of solutions (bestx6(:,2)). It will be required for me for further research and optimisation.

As i followed your steps .. it will be easy if you can suggest the rest of it.

Thanks in advance.

[Walter Roberson]

full_sol = subs(sol, {x4, x6}, {x4_values, X6(:,2)});

The result will be a struct with arrays x1, x2, x3, x5, each of which has one entry for each time value.

You might need to work on that a bit to get the order right (that is, you might need to substitute x6 before x4.)

for K = 1 : length(x4_values)
    x4 = x4_values(K);
end 

That loop is useless. It just gives you x4 = x4_values(end) .

[Jhon Paul]

Thank you veyr much .. I have got the whole thing ...

[Jhon Paul]

I have got the solutions as you have suggested for T = 2000:500:18000;

x0 = [603, 5*(10^4), 0.011, 0, 1003, 1373];

x4_values = x0(4) + 2.2*10^-6 * T;

for K = 1 : length(x4_values)

x4 = x4_values(K);

end

x4= double(x4_values)

x61= real(vpa(subs(bestx6(1),x4)));

Then I followed to solve each variable from the equations:

x2=(1373-x61)./(0.12*x4)

x3 =(2.5*10E-4*x61.^2)./x2

x1 = (70600+x2)./200

However in the actual problem, there is variable velocity term ( for 2.2*10^-6) which needs a loop. I have created the pseudocode, if you can help. But I am stuck in indexing the variables for the loop.

T = 2000:500:18000;
x0 = [603, 5*(10^4), 0.011, 0, 1003, 1373];
  x1 = zeros(1,13);
  x2 = zeros(1,13);
  x3 = zeros(1,13);
  x4 = zeros(1,13);
  x5 = zeros(1,13);
  x6 = zeros(1,13);
for K = 1 : length(T)
      v(k) = ((147408*x6.^2).*exp(-39211/x6))./x2;
      x4_values = x0(4) + v(k+1)* T;
      x61= real(vpa(subs(bestx6(1),x4_values)));
      x2=(1373-x61)./(0.12*x4);
      x3 =(2.5*10E-4*x61.^2)./x2;
      x1 = (70600+x2)./200
end