MATLAB Answers

Fsolve with Loop and to store variable

3 views (last 30 days)
Jhon Paul
Jhon Paul on 22 Jan 2021
Commented: Jhon Paul on 31 Jan 2021
I am having the following set of 4 non linrear equations with 6 variable x1,x2,x3,x4,x5,x6
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;
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 !

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 = [
2*x2*(15*10^-3)-x3-x4-x5+x1==0 %equation was corrupt, had to guess what it was
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
Jhon Paul on 31 Jan 2021
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);
x4= double(x4_values)
x61= real(vpa(subs(bestx6(1),x4)));
Then I followed to solve each variable from the equations:
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)));
x3 =(2.5*10E-4*x61.^2)./x2;
x1 = (70600+x2)./200

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!

Translated by