Why is there a significant difference in the assignment results of functions using "subs" and "feval"

I need a large amount of data for assignment, and I need to use the handle function to accelerate it, but the result calculated using "feval" is unacceptable

clear,clc
syms z
load FF.mat Th3
z3 = 0.0016:0.00022:0.006;
T1 = vpa(subs(Th3,z,z3));
disp(vpa(T1.',14))

Th3f = matlabFunction(Th3);
T2 = feval(Th3f,z3);
disp(T2.')

    3.5028
    3.9585
    4.4094
    4.8132
    5.1376
    5.3597
    5.4692
    5.4640
    5.3344
    5.0730
    4.6994
    3.9508
    2.9207
    0.9343
   -2.8567
   -7.1941
  -31.3176
  -55.4727
 -173.9580
 -324.4082
 -787.7051

[Dyuman Joshi]

"but the result calculated using "feval" is unacceptable"

Why/How exactly is it unacceptable?

[sen]

According to the formula for solving, the last result needs to be close to 0, while the calculated result for “feval” is -787.7051

[sen]

Thank you for your reply. Could you please run the code I submitted? I would like to know how to use "feval" to obtain results that are close to those obtained by "subs"

Answer 1

[Walter Roberson]

syms z
load FF.mat Th3
z3 = 0.0016:0.00022:0.006;
T1 = vpa(subs(Th3,z,z3));
disp(vpa(T1.',14))

Th3f = matlabFunction(Th3);
T2 = feval(Th3f,z3);
format long g
disp(T2.')

          3.50282686149745
          3.95854248851538
          4.40936088562012
           4.8132483959198
          5.13755702972412
          5.35973119735718
          5.46921586990356
          5.46400547027588
          5.33437728881836
          5.07296371459961
          4.69937896728516
          3.95075988769531
          2.92066955566406
            0.934326171875
           -2.856689453125
           -7.194091796875
          -31.317626953125
              -55.47265625
           -173.9580078125
            -324.408203125
            -787.705078125

T12 = double(T1(:)) - T2(:);
disp(T12)

      7.52585900190006e-05
      0.000255996608482079
      0.000548036419161235
       0.00103877265106078
       0.00206573817902189
       0.00575961188859786
        0.0133017952813956
        0.0256612421124434
        0.0564997901084956
         0.119048502399311
         0.199439386434167
         0.565610566987848
          1.12954361764356
          2.57470185401177
          5.76486565117093
          9.46704282663683
           32.957595256356
          56.5276612838938
          174.523944755596
          324.618623611357
             787.705078125

plot(z3, T12)

[Walter Roberson]

So the difference in values gets worse for larger z3 values.

Your expression is complicated enough that trying to simplify it takes more than 4 minutes.

It is very typical that complicated expressions lose precision when evaluated numerically. The only way to get around that is to rewrite the expressions to try to maximize precision over the range of interest.

[sen]

This is already the simplest form I can do. I have only expanded the first term of the series. Is there any other way to accelerate the "subs"

[sen]

And what are the reasons that lead to different results for the "feval"

[Walter Roberson]

syms z
syms delta real
load FF.mat Th3
z3 = 0.0016:0.00022:0.006;
T1 = subs(Th3, z,z3(1));
T1d = subs(Th3, z, z3(1)+delta);

vpa(T1, 10)

ans = 

3.502901928

vpa(T1d, 10)

ans = 

The calculation involves sums of negative and positive exponentials with a difference of δ in the input value making a difference of more than

Order of calculations would obviously matter a fair bit.

[sen]

Even very small fluctuations can cause significant differences in results; The calculation time of "subs" is very long, and there are problems with the result of the calculation of "feval"; How can I improve the algorithm?

I want to assign values quickly and accurately, what should I do

[Walter Roberson]

syms z delta real
load FF.mat Th3
Th3f = matlabFunction(Th3);
format long g

T3 = feval(Th3f, z);

T3D = vpa(Th3 - T3, 14)

T3D = 

subs(T3D, z, 0.0016)

ans = 

0.000069891907750774627725761978241548

subs(T3D, z, 0.006)

ans = 

700.41596617851346360052434728079

What this tells you is that after the symbolic expression is converted to numeric function, if you evaluate the numeric function at a symbolic parameter to recover a symbolic expression, and you take the difference between that and the original symbolic expression... then the difference is pretty significant near the upper end of your range.

You are losing a notable amount of precision in converting to numeric form.

[Walter Roberson]

To meet your goals, you are probably going to have to rent time on a supercomputer.

Your expressions are highly non-linear.

[sen]

Thank you very much for your help, it is indeed so;

Can only use the "subs"？it's too slow

[Walter Roberson]

fprintf('step 1\n');

step 1

syms z delta real
fprintf('step 2\n');

step 2

load FF.mat Th3
fprintf('step 3\n');

step 3

ch = children(Th3);
fprintf('step 4\n');

step 4

ch = [ch{:}].';
fprintf('step 5\n');

step 5

zch = subs(ch, z, 0.006);
fprintf('step 6\n');

step 6

[~, order] = sort(double(zch));
fprintf('step 7\n');

step 7

parts = vpa(zch(order), 14)

parts = 

sum(parts)

ans = 

You are adding together a number of components in the range +/- 1e17 and expecting to get a result on the order of 1e-6 -- which is basically noise level compared to the range of values you are dealing with.

[sen]

Thank you for your reply. I understand what you mean.

Actually, many times the boundary conditions I define are very small, such as 0; However, during the process, a large number of values appeared, which should be related to the material parameters I calculated, with some smaller values being used as denominators;

And this is not a separate equation. I have a set of similar equations that can lead to this situation. I gave a simple equation as an example.

At the same time, sometimes I need to perform FFT operations on the calculated results, and such a large loss of accuracy can lead to problems.

I don't want to add some data correction in the program because it is difficult to apply to various types of situations. Therefore, I want to start from the program and accelerate the assignment that can preserve accuracy.

I want to know if there are other methods besides the "subs" that can calculate faster and reduce my accuracy loss

I just want to calculate faster

[Walter Roberson]

"I want to know if there are other methods besides the "subs" that can calculate faster and reduce my accuracy loss"

No. Your values range between +/- 10^17 and you want the total to be accurate to better than 10^-6. In order to have a hope of that, you need your calculations to be accurate to at least 25 digits, which is something that is not possible using double precision numbers.

A few weeks ago I did see someone post reference to a replacement symbolic toolbox that they claimed was much faster. It will take me some time to locate that discussion again.

[Walter Roberson]

Found it: https://www.mathworks.com/matlabcentral/answers/382123-vpa-4503599627370491-5#answer_1304861

[sen]

Thank you very much for your help