MATLAB Answers

Fsolve with Loop and to store variable

65 views (last 30 days)
Jhon Paul
Jhon Paul on 22 Jan 2021
Commented: Jhon Paul on 31 Jan 2021 at 18:24
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 !

  3 Comments

Walter Roberson
Walter Roberson on 22 Jan 2021
2*x2*(15*10^-3)-x3-x4-x5=x1=0
please recheck that. It implies that x1 must always be 0.
Walter Roberson
Walter Roberson on 22 Jan 2021
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.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 22 Jan 2021
Edited: Walter Roberson on 22 Jan 2021
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.

  20 Comments

Walter Roberson
Walter Roberson on 29 Jan 2021 at 22:24
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
Jhon Paul on 30 Jan 2021 at 9:45
Thank you veyr much .. I have got the whole thing ...
Jhon Paul
Jhon Paul on 31 Jan 2021 at 18:24
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

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!