How can I correctly substitute symbolic variables with decimal numbers like 2.543342? The subs function can't do it properly...

This is the function I wrote:

function [ nsquare ] = sellmeier_formula( type , coefficients )
%SELLMEIER_FORMULA This functions gives the sellmeier formula as a symbolic
%function of lambda in micrometers.
%   Inputs:
%       
%   type - It's the type of sellmeier formula to be used.
%   
%   coefficients - It's a vector containing the coefficients to be used for
%   the sellmeier formula.
% 
%   Outputs:
% 
%   nsquare - It's a symbolic function with argument lambda of the square
%   of the refractive index n corresponding to the type of sellmeier
%   formula used.
%   !!!WARNING!!!LAMBDA NEEDS TO BE EXPRESSED IN UNITS OF MICROMETERS
%       

%Input arguments check
minargs=2;
maxargs=2;
narginchk(minargs, maxargs)

%Input data class and validity check
if ~isnumeric(type)
    warndlg('type must be a numeric input');
    error(['type has to be a numeric input']);
elseif ~isnumeric(coefficients)
    warndlg('coefficients must be a numeric input');
    error(['coefficients has to be a numeric input']);
end

% Function beginning:

lambda = sym('lambda');
D = zeros(1,17);
C = sym('c',[1 17]);

for i=1:length(coefficients)
    D(i) = coefficients(i);
end

for i=1:17
C(i) = subs(C(i),D(i))
end

switch type
    case 1
        nsquare = 1+C(1)+(C(2).*lambda.^2)./(lambda.^2-C(3).^2 )+(C(4).*lambda.^2)./(lambda.^2-C(5).^2 )+(C(6).*lambda.^2)./(lambda.^2-C(7).^2 )+(C(8).*lambda.^2)./(lambda.^2-C(9).^2 )+(C(10).*lambda.^2)./(lambda.^2-C(11).^2 )+(C(12).*lambda.^2)./(lambda.^2-C(13).^2 )+(C(14).*lambda.^2)./(lambda.^2-C(15).^2 )+(C(16).*lambda.^2)./(lambda.^2-C(17).^2 );
    case 2
        nsquare = 1+C(1)+(C(2).*lambda.^2)./(lambda.^2-C(3) )+(C(4).*lambda.^2)./(lambda.^2-C(5) )+(C(6).*lambda.^2)./(lambda.^2-C(7) )+(C(8).*lambda.^2)./(lambda.^2-C(9) )+(C(10).*lambda.^2)./(lambda.^2-C(11) )+(C(12).*lambda.^2)./(lambda.^2-C(13) )+(C(14).*lambda.^2)./(lambda.^2-C(15) )+(C(16).*lambda.^2)./(lambda.^2-C(17) );
    case 3
        nsquare = C(1)+C(2).*lambda.^(C(3) )+C(4).*lambda.^(C(5) )+C(6).*lambda.^(C(7) )+C(8).*lambda.^(C(9) )+C(10).*lambda.^(C(11) )+C(12).*lambda.^(C(13) )+C(14).*lambda.^(C(15) )+C(16).*lambda.^(C(17) );
    case 4
        nsquare = C(1)+(C(2).*lambda.^(C(3) ))./(lambda.^2-C(4).^(C(5) ) )+(C(6).*lambda.^(C(7) ))./(lambda.^2-C(8).^(C(9) ) )+C(10).*lambda.^(C(11) )+C(12).*lambda.^(C(13) )+C(14).*lambda.^(C(15) )+C(16).*lambda.^(C(17) );
    case 5
        nsquare = C(1)+C(2).*lambda.^(C(3) )+C(4).*lambda.^(C(5) )+C(6).*lambda.^(C(7) )+C(8).*lambda.^(C(9) )+C(10).*lambda.^(C(11) );
    case 6
        nsquare = 1+C(1)+C(2)./(C(3)-lambda.^(-2) )+C(4)./(C(5)-lambda.^(-2) )+C(6)./(C(7)-lambda.^(-2) )+C(8)./(C(9)-lambda.^(-2) )+C(10)./(C(11)-lambda.^(-2) );
    case 7
        nsquare = C(1)+C(2)./(lambda.^2-0.028)+C(3).*(1./(lambda.^2-0.028)).^2+C(4).*lambda.^2+C(5).*lambda.^4+C(6).*lambda.^6;
    case 8
        nsquare = 3./(1-(C(1)+(C(2).*lambda.^2)./(lambda.^2-C(3) )+C(4).*lambda.^2))-2;
    case 9
        nsquare = C(1)+C(2)./(lambda.^2-C(3) )+(C(4).*(lambda-C(5)))./((lambda-C(5)).^2+C(6) );
    otherwise
        nsquare = -1;
end
if D(1)==0
    D = D(2:end);
    D = D(D~=0);
elseif D(1)~=0
    D = D(D~=0);
end

end

But when using the subs function in the for loop I get fractional numbers that are not correct or huge numbers that have nothing to do with the original decimal number... Can anybody help with this problem please?

Answer 1

[John D'Errico]

Sure it does work properly. It is just that a decimal number is not exactly what you think it is. Floating point numbers are not represented exactly in the binary form used to store a number.

x = sym(2.543342)
x =
2863548520868933/1125899906842624

So MATLAB sym converts it to a fraction, composed of the ratio of integers. Those huge numbers have EVERYTHING to do with the original number. If you want to see the value stored, then use vpa.

vpa(x)
ans =
2.5433419999999999916440174274612

If you REALLY want to store the exact decimal value, then do it like this:

x = sym('2.543342')
x =
2.543342

subs is no different. It uses the same tools.

[test run]

Hello John D'Errico , I haven't tried the

vpa(x)

option yet , but I have tried allready the

x = sym('2.543342')

The problem using this is that later I will have the need for those coefficients to be actual numbers of class double in order to use the expression of nsquare for different values of lambda. You meant to say that by leaving the coefficients that way

for i=1:17
C(i) = subs(C(i),num2str(D(i)));
end

then I'll be able to do numeric calculations when I replace lambda by it's numerical value???

thanks for your time in advance

[John D'Errico]

NO. I did not say that. You need to understand that the ratio of integers IS the same as the number that you substituted. If you want to convert a result to look like it is in decimal form, then you can use vpa. Or you can use double to convert the result to double precision.

What I am trying to tell you is that

2.543342

is NOT in fact stored as 2.543342, when it is used as a number. In fact, you cannot store 2.543342 as a decimal number in double precision, because binary storage is used.

x = sym(2.543342);
vpa(x,100)
ans =
2.54334199999999999164401742746122181415557861328125

I've shown you the closest binary number to the decimal number you want to use. In fact, when you type 2.543342 in matlab, you get the number seen above.

So the symbolic toolbox stores that as a ratio of integers. It can work with integers very well. But floating point numbers are difficult, since they always lose some precision in the least significant bits.

So you can use subs just as you did before. There is no need to convert a double into a string to use subs.

For example:

syms x
y = x + 2;
z = subs(y,x,3.14159)
z =
5788915702022967/1125899906842624

MATLAB stores the result as a ratio of integers. No problem converting it into a double.

double(z)
ans =
                   5.14159

But is that truly 5.14159?

vpa(z,100)
ans =
5.14158999999999988261834005243144929409027099609375

No. In fact, it is not exactly what you think it is, but an approximation that looks like you expect, if you rounded it after about 16 decimal digits or so.

In fact, even floating point numbers in symbolic form have limits.

z = subs(y,x,'3.14159')
z =
5.14159

vpa(z,100)
ans =
5.141589999999999999999999999999999999999838251977326853238903501852478218994732151619388160757040396

So passing it in as a string was still represented using an approximation. (I can't recall the number of binary bits used here. Easy to figure out though.)

To get the full true precision, I could have done this:

z = subs(y,x,'314159/100000')
z =
514159/100000

vpa(z,100)
ans =
5.14159

The point is, a ratio of integers is not a problem as long as you understand to use double or vpa. Regardless, you need to understand how numbers are stored.

[test run]

Thank you for your patience and explanations I understand now what you explained to me. From your answer I see what I have to do and also using vpa helps showing the numbers with a decimal representation. But also the floating point number approximation as a fraction of 2 integers gives a very accurate result. Thank you again for the explanation.

Answer 2

[Mischa Kim]

Are you referring to the C(i)?

How about

C(i) = vpa(subs(C(i),D(i)))

[test run]

Thank you that helped me show the numbers in a decimal representation.

Answer 3

[Walter Roberson]

http://www.mathworks.com/help/symbolic/sym.html#input_argument_flag

'r' (default) | 'd' | 'e' | 'f'

"When sym uses the rational mode, it converts floating-point numbers obtained by evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q, and 10^q for modest sized integers p and q to the corresponding symbolic form. This effectively compensates for the round-off error involved in the original evaluation, but might not represent the floating-point value precisely. If sym cannot find simple rational approximation, then it uses the same technique as it would use with the flag 'f'."

If you want the symbolic number to be exactly the same value as the binary number you are substituting, then you need to use the 'f' flag:

"When sym uses the floating-point mode, it represents all values in the form N*2^e or -N*2^e, where N >= 0 and e are integers. For example, sym(1/10,'f') returns 3602879701896397/36028797018963968 . The returned rational value is the exact value of the floating-point number that you convert to a symbolic number."

For example,

x = 2.543342;      %x becomes a binary floating point number, not a decimal number!
xr = sym(x);       %same as sym(x,'r'), xr becomes a symbolic integer ratio that _approximates_ the binary floating point number
xf = sym(x, 'f');  %xf becomes a symbolic integer ratio whose value is exactly the same as the binary floating point number

[test run]

Thank you for the explanations it helped me understand.