How to plot graph of value of variable when a parameter is varied

Hi all,

I have a system of equations and I would like to solve for e, τ, and k. However, I want to assume that the parameter R can take on several different values (e.g. 0.5, 1, 1.5). In particular, I hope to plot a graph of how the solution for any of my variables e, τ or k varies with R. I am really struggling with how this can be implemented in MATLAB. Any help will be much appreciated.

Answer 1

[Ameer Hamza]

See fsolve(): https://www.mathworks.com/help/optim/ug/fsolve.html

Write your system of equation in the form

Where , and are the three equations in your queston. Then apply fsolve() function.

This is the code:

f_1 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - 0.336*e.^-0.3./(tau+R.*k);
f_2 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - 0.5./tau;
f_3 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - (R-0.3)./(tau+R.*k);

f = @(e, tau, k, R) [f_1(e,tau,k,R); f_2(e,tau,k,R); f_3(e,tau,k,R)];

R = [0.5, 1, 1.5];
sol = zeros(numel(R), 3);
for i=1:numel(R)
    sol(i,:) = fsolve(@(x) f(x(1),x(2),x(3), R(i)), rand(1,3));
end

plot(R, sol);
legend({'e', 'tau', 'k'});

[Teck Lim]

Thanks for your quick reply. It is immensely helpful.

I tried out the code in MATLAB and have some questions:

1) Am I right to say that the linear plots for e, tau and k I see are just extrapolations from the solutions at the 3 points R=0.5,1,1.5? Meaning to say MATLAB does not actually calculate the value of the solution at all points in the domain from 0.5 up to 1.5.

2) For some reason the solutions that the algorithm generates is much higher than what I got from pen and paper calculation. All my manual solutions are below 1 but the algorithm gives me solutions from 45 all the way up to 85. Am I able to impose extra conditions such as having the algorithm check for optimal solutions at smaller values?

[Ameer Hamza]

1) Yes, MATLAB does not automatically interpolate the values. These are just straight lines connecting 3 points.

2) Does this system of equations have multiple solutions? Can you share one of the solution you got from pen and paper calculations?

[Teck Lim]

For my pen and paper calculations, I assumed R=1 and I obtained:

τ=0.01698

e=0.0866

k=0.00679

Certainly, I am not ruling out multiple solutions. How would the algorithm respond if that was the case? Does that mean the solution will depend on the initial guess?

[Ameer Hamza]

Teck, your equation does seem to have multiple solutions, and the solution you mentioned for R=1 is correct. I tried different initial points, and values of R. The fsolve() is able to find a solution if R >= 0.8, but for smaller values, it cannot find a solution

f_1 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - 0.336*e.^-0.3./(tau+R.*k);
f_2 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - 0.5./tau;
f_3 = @(e, tau, k, R) 1./(0.8*e.^0.7-e-tau-k) - (R-0.3)./(tau+R.*k);

f = @(e, tau, k, R) [f_1(e,tau,k,R); f_2(e,tau,k,R); f_3(e,tau,k,R)];

R = [0.5, 0.8, 1, 1.5, 2];
sol = zeros(numel(R), 3);
for i=1:numel(R)
    sol(i,:) = fsolve(@(x) f(x(1),x(2),x(3), R(i)), 1e-6.*rand(1,3));
    f(sol(i,1), sol(i,2), sol(i,3), R(i)) % to check if the 3 equations are true.
end
plot(R, sol);
legend({'e', 'tau', 'k'});

[Teck Lim]

Thanks again for your time.

Answer 2

[Torsten]

Your equations can only be solved for one value of R, namely R = 0.3 + 0.366*exp(-0.3).

This can be seen by dividing equation 1 by equation 3.

[Ameer Hamza]

Torsten, e is not Euler's number in this question. OP used e to denote a variable.

[Torsten]

ok, then at least we can explicitly solve for e depending on R.

Taking the reciprocals of all the 6 expressions and inserting the result for e from the first step, we can solve for tau and k since we have a linear system of equations.

[Teck Lim]

Thanks for the suggestion. I will try that out.

[Ameer Hamza]

Torsten's analysis is correct. This system does have a closed-form solution. I used a symbolic toolbox to solve this system of equation, and it gives more accurate results (as expected) at R=0.5 as compared to the numeric solver. Following code shows the solution with symbolic approach

syms e tau k R

eq1 = 1./(0.8*e.^0.7-e-tau-k) - 0.336*e.^-0.3./(tau+R.*k);
eq2 = 1./(0.8*e.^0.7-e-tau-k) - 0.5./tau;
eq3 = 1./(0.8*e.^0.7-e-tau-k) - (R-0.3)./(tau+R.*k);

r = [0.5 1 1.5];
sols = zeros(numel(r), 3);
for i=1:numel(r)
    sol = solve([eq1==0 eq2==0 eq3==0], [e tau k]);
    sols(i,1) = subs(sol.e(1), R, r(i));
    sols(i,2) = subs(sol.tau(1), R, r(i));
    sols(i,3) = subs(sol.k(1), R, r(i));
end

sols = double(sols);
plot(r, sols);
legend({'e', 'tau', 'k'});