Can't find error in code.

I've thrown together something from various forums in an attempt to solve a numerical problem. I'm trying to find the values for N and V for a range of n. There's a syntax error on line 3 and a parse error in line 5, possibly due to incorrect parentheses. I've checked the parentheses many times and I believe they're fine, so I'm not sure what the issue is. Here's the code:

syms N V

for n = 1:5
  [solN,solV] = solve(N == 51.4 + log(V^(1/4)*(2*n*(2*n-1))^(n/4)/(10^(16)))
+(1/2)*log((2*N+2*n-1)/(2*n-1)), 
  V == ((5.9*10^(33)*n)/(2*n*(2*N+2*n-1))^((n+1)/2))^2)
  disp(n)
end

Any help would be appreciated.

Answer 1

[David Goodmanson]

Hi David, [This is a new answer since modifying my original erroneous answer may not be the most efficient way to go.]

The equations do appear to have good numerical solutions.

Ns = zeros(5,1);
Vs = zeros(5,1);
errel = zeros(5,2);        % relative error in equations for N and V
g = @(N,n) ((5.9e33*n)./(2*n*(2*N+2*n-1)).^((n+1)/2)) ...
    -(exp(2*(N-51.4))./(((2*N+2*n-1)/(2*n-1))*((2*n*(2*n-1))^(n/4)/1e16)^2));
for n = 1:5
  N = fzero(@(N) g(N,n),50);
  Ns(n) = N;
  V = ((5.9e33*n)./(2*n*(2*N+2*n-1)).^((n+1)/2))^2;
  Vs(n) = V;
  errel(n,1) = 1 - ( 51.4 + log(V^(1/4)*(2*n*(2*n-1))^(n/4)/(10^(16))) ...
    +(1/2)*log((2*N+2*n-1)/(2*n-1)) )/N;
  errel(n,2) = 1 -  ( ((5.9*10^(33)*n)/(2*n*(2*N+2*n-1))^((n+1)/2))^2 )/V;
end
format long
Ns
Vs
errel
Ns =
 53.2655
 52.2693
 51.6005
 51.0656
 50.6132

Fzero is looking at the difference between two expressions for sqrt(V). Vs does not display very well but of course the individual elements are available.

[David Anthony]

Excellent!

Thanks a lot for your time, those values look very promising. Now to teach myself what you wrote actually means.

[David Anthony]

Hi again. Upon closer inspection I noticed some errors I made with the original expressions. More generally, as a beginner I don't understand what your code is doing (especially what the g is doing), so if possible I'd prefer something simpler that I can generalise. I've kept your answer as 'solved' because it is solved, and the further issues are my own. All that said, can I ask why the following doesn't work?

syms N V
for n = 1:5
    [solN,solV] = solve(N == 50.37+log(V.^(0.25)*(2*n(2*n-1).^(n./4))./(10^(16)))...
        +(n./2)*log((2*N+2*n-1)./(2*n-1)),...
        V == (3.495*10^(67)*n.^2)./((2*n(2*N+2*n-1))^(n+1)))
end

I get the following error:

Error using sym/subsindex (line 663)
Invalid indexing or function definition. When defining a function, ensure that the body of the function is a SYM object. When indexing, the input must be numeric, logical or ':'.

Error in chaotic3 (line 3)
    [solN,solV] = solve(N == 50.37+log(V.^(0.25)*(2*n(2*n-1).^(n./4))./(10^(16)))...

Thanks again for the help so far, and any more if possible would be very much appreciated.

[David Goodmanson]

Hi David,

Looks like the long format versions of N didn't make it into my post after all but now it does not matter.

I can't say anything about the syms solution (except maybe it wants n to be a symbolic variable too) because unfortunately I don't have any toolboxes. But because of Matlab for home use which is my situation I will now be able to purchase some.

I can say something about the solution above, though. It would have been much clearer to define g in an mfile, because then you can abbreviate expressions to make it easier to solve the equations algebraically. I did that, but then for whatever reason put it back inline with everything expanded out so it's not so easy to read. So as it stands, it's not good Matlab code. Here is an mfile with some annotations:

function D = diffsV(N,n)
% difference between expressions for sqrt(V) from two different functions.
% both are satisfied when D = 0.
%  
% N == 51.4 + log(V^(1/4)*(2*n*(2*n-1))^(n/4)/(10^(16))) ...
%      +(1/2)*log((2*N+2*n-1)/(2*n-1))
% V == ((5.9*10^(33)*n)/(2*n*(2*N+2*n-1))^((n+1)/2))^2

C = 51.4;
E = 1e16;
B = 5.9e33;
a = (2*n*(2*n-1))^(n/4);
FN = (2*N+2*n-1)/(2*n-1);
GN = (2*n*(2*N+2*n-1))^((n+1)/2);

% V = ((B*n)/GN)^2
sqrtV1 = (B*n)./GN;

% N = C +      log(V^(1/4)*(a/E)  ) + (1/2)*log(FN)
% N = C + (1/2)log(V^(1/2)*(a/E)^2) + (1/2)*log(FN)
% invert this
sqrtV2 = exp(2*(N-C))./(FN.*(a/E).^2);

D = sqrtV1 - sqrtV2;

It's certainly not totally generalizable, but I hope that helps.

I solved the second equation for sqrt(V) not by symbolic toolbox but with pencil and paper. We used to do that a lot when I was in school. (Seriously, symbolic computation is a great advance, but while a great deal was gained, something was lost too. Overreliance on syms is a problem.)

[David Anthony]

Thanks, using that breakdown of the equations I can get my code to work. No idea what the issue was, but doesn't matter now.

Answer 2

[Star Strider]

When in doubt with a long expression, break it up into its component expressions and check each for matching prens.

This works. I will leave it to you to determine if it gives the correct result:

syms N V

for n = 1:5
    Term1 = log(V^(1/4)*(2*n*(2*n-1))^(n/4)/(10^16));
    Term2 = (1/2)*log((2*N+2*n-1)/(2*n-1));
    solN = solve(N == 51.4 + Term1 + Term2, N);

      V = ((5.9*10^(33)*n)/(2*n*(2*N+2*n-1))^((n+1)/2))^2;
      Vv(n,:) = subs(V, N, solN);
      disp(n)
  end
  Vv = vpa(Vv, 5)

The loop is causing problems with the code as your originally wrote it, so my changes were necessary for it to run in the loop.

[David Anthony]

Thank you.

Can I ask what Vv is doing? I'm trying to solve for V and N simultaneously for the given range of n - is that what the code is doing at present? If I put a disp(solN) disp(solV) in the loop, I get values of N that are expected for n=2:5, but none N for n=1, and also no displayed values for any V. I do get four of the five values roughly in line with what I'd expect for V after the loop, though.

[Star Strider]

My pleasure.

When I corrected my code to understand that you wanted to continue the first line, I ran the loop and kept getting:

Warning: Cannot find explicit solution. 
> In solve (line 316)

Since this is not uncommon, (and since I don’t know what problem you’re solving, I didn’t want to use 'IgnoreAnalyticConstraints',true). I decided to give it to fsolve instead.

The Code —

VN = @(N,V,n) [51.4 + log(V^(1/4)*(2*n*(2*n-1))^(n/4)/(10^(16))) +(1/2)*log((2*N+2*n-1)/(2*n-1)), ...
    ((5.9*10^(33)*n)/(2*n*(2*N+2*n-1))^((n+1)/2))^2];
for n = 1:5
    VNsol = fsolve(@(b) VN(b(1),b(2),n), rand(2,1));
    Nv(n,:) = VNsol(1);
    Vv(n,:) = VNsol(2);
    disp(n)
end

fprintf(1,'\t\tN\t\tV\n')
fprintf(1,'\t%7.2f\t%7.4f\n', [Nv, Vv]')

      N      V
   154.30   0.6335
    67.57   0.5647
    35.15   0.1433
    28.87   0.8367
    26.14   0.3210

This starts with a new, arbitrary initial parameter estimate vector each time, and the results — while approximately the same — never are the same, so you may want to experiment with different initial estimates. All I can verify is that the code works. I don’t know if it produces the results you want. (Since I actually have no idea what problem you are solving, I can’t apply what expertise I may have to it.)

[David Anthony]

Thanks again.

I've copied your code and it doesn't seem to work on my end. It says:

"Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead. > In fsolve at 285 In chaotic2 at 4

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient, but the vector of function values is not near zero as measured by the default value of the function tolerance."

That's without with fprint. With the fprint I get the error:

"Error using fprintf Function is not defined for 'sym' inputs.

Error in chaotic2 (line 11) fprintf(1,'\t%7.2f\t%7.4f\n', [Nv, Vv]')"

For the record the values you've shared aren't looking promising for me. The N values are supposed to be less than 60ish, and the V values are very roughly supposed to be of order 10^60... That's for me to worry about, though!

[Star Strider]

My pleasure.

I am using R2016b. There could be version differences.

Run my code with the symbolic code deleted or commented-out. The code I posted worked for me (as I documented). I don’t know the problem you’re solving, so I have no suggestions on how best to solve it.

The initial parameter estimates are important. See if:

VNsol = fsolve(@(b) VN(b(1),b(2),n), [50, 1E+60]);

and:

fprintf(1,'\t%7.2f\t%7.4E\n', [Nv, Vv]')

work.

That is the best I can do.

[David Anthony]

Thanks a lot for your time.

[Star Strider]

As always, my pleasure.

Answer 3

[David Goodmanson]

Hi David,

Looks like the parentheses for line 3 aren't correct after all. Needs one more parenthesis at the end of the statement. And V has one too many, not needing the parenthesis at the end..

[David Anthony]

Thanks for the reply -

Currently the parenthesis at the end of line 5 is closing the 'solve' on line 3. If I add an extra parenthesis at the end of the statement, the 'solve' ends before the V term. Regardless, if I make the change (unless I've misunderstood you) I have the same error. Apologies if I'm misunderstanding you.

[David Goodmanson]

I didn't have the right idea and I see what you mean. At this point could you stand to use the three-dot continuation ... at the end of line 3?

[David Anthony]

Ah, I didn't realise that new paragraphs actually impacted the code in any way, I added them just for ease of reading.

Putting everything on one line fixes the errors, but results in an emptysym for N, and no explicit solutions for V. This may be because the equations don't have an analytical solution, but I'm not sure how to fix it.