Unable to plot a complex function with EZplot

I am trying to plot this function (fun):

with r as the only variable over a set of defined constants.

I have already produced a plot of this with a given set of constants:

a=1.59e-04;
C0=75.5;
C1=80.4;
D=5.7e-07;
t=1995840;
syms n;
syms r;
fun = @(r) C1+(C0-C1).*(1+((2.*a)./(r.*pi)).*vpa(symsum((((-1).^(n))./n).*sin((n.*r.*pi)./a).*exp(((-D.*n.^(2).*pi.^(2).*t)./a.^(2))),n,1,inf)));
ezplot(fun, [0,a])

The plot appears correct, nevertheless several errors come up:

Error using message
In 'MATLAB:axis:UnknownOption', data type supplied is incorrect for parameter {1}.
Error in axis (line 176)
              error(message('MATLAB:axis:UnknownOption', cur_arg));
Error in ezplot1ezplot1 (line 548)
      axis(cax, [xmin xmax ymin ymax]);
Error in ezplot (line 144)
              [hp, cax] = ezplot1(cax, f{1}, vars, labels, args{:});

After changing these constants:

a=8.06e-04;
C0=42.35;
C1=46;
D=1.1e-17;
t=864000;
syms n;
syms r;
fun = @(r) C1+(C0-C1).*(1+((2.*a)./(r.*pi)).*vpa(symsum((((-1).^(n))./n).*sin((n.*r.*pi)./a).*exp(((-D.*n.^(2).*pi.^(2).*t)./a.^(2))),n,1,inf)));
ezplot(fun, [0,a])

the program will run for an extended period of time (>1 day) and no plot is created. I assume there may be an issue in symbolizing n or r (i.e. operation running on non-real variables).

Answer 1

[John BG]

James

the point that you are plotting with ezplot means you are using a discrete time base.

Try truncating the sum series

where whether 1/n is small enough that the terms do not add anything significant or exp(-n) also reducing the terms to zero.

Would it be possible to tell what exactly, what variables are you plotting in the plot of the question?

John

[e-mail]

[James Valentine]

Hello John,

The function refers to 3d diffusion of a sphere of radius a of initial concentration C 1 in an infinite medium of surface concentration Co. The plot shows distance along the diameter of the sphere (r) against equilibrated concentration (C) at a discrete point in time (as you were mentioning). Change in C with respect to r is dependent on D, the diffusion coefficient (in units of length^2/time, dependent on the type of material involved).

I realize truncating the series would have a negligible effect on the values, but I was aiming for plotting with respect to the infinite sum. Nonetheless, this should be feasible (as 1/n will become increasingly small regardless).

James

Answer 2

[Walter Roberson]

My investigations suggest that plot is not correct.

Look at the exp(-n^2*Pi^2*D*t/a) term. Let P = D*t/a. If P > 1, then exp(-n^2*Pi^2*P) gets small quickly as n increases. Experimentally, when P >= 4, the rate at which it gets small is sufficiently high that the sum becomes entirely dominated by the first term, when n = 1, forming the expression

C1+(C0-C1)*(1-2*a*sin(r*Pi/a)*exp(-D*Pi^2*t/a^2)/(r*Pi))

Your D*t/a is 2069760000000000000/45995321621 which is far greater than 4. You effectively have about exp(-4.4E8) as your multiplier of sin(r*Pi/a), which is about 0.5525277122674722e-192881416 . As close to 0 as almost makes no difference. But what is it being multiplied by? By 2*a*sin(r*Pi/a) . Your a is far less than 1/2. sin(r*Pi/a), whatever it is, is going to be between -1 and +1 . So with your exp() term being effectively 0, your 2*a*sin(r*Pi/a)*exp(-D*Pi^2*t/a^2)/(r*Pi) is going to be effectively 0. And that's going to leave C1+(C0-C1) = C0 as your result for numeric purposes -- for all values of r. Your plot is going to be a constant plot.

Your D or your t or your a would have to be greatly different in order to change this result.

[James Valentine]

Hello Walter,

Both t and D (for my purposes) are plotted on the order of milliseconds (with a being micron-scale) with respect to time, so I understand that subtle variations in these variables will be equivalent to plotting nearly identical solutions.

James

[John BG]

the expression inside the sum series:

(-1).^n./n*sin(r/a*pi*n).*exp(-D/a^2*pi^2*n^2*t)=
(-1).^n./n*sin(r/a*pi*n).exp(-D/a^2*pi^2*n^2)^t

decreases quickly with n.

Please correct if wrong but n=1:1:Inf n does not take values 0 or negative, precisely because you said it's a ball. So

format long
n=1:10
exp(-n.^2)
=

     0.367879441171442   0.018315638888734   0.000123409804087   0.000000112535175

     0.000000000013888   0.000000000000000   0.000000000000000   0.000000000000000

     0.000000000000000   0.000000000000000

no matter how big exp(t) is that once multiplied by this coefficients it all goes to zero.

You mention it is a spherical diffusion that in all applications I have read or seen: wireless radiation, heat, sound, chemical explosion, it all quickly decreases with 1/R 1/R^2 1/R^3 (electromagnetic radiated far field or near field) unless you have some kind of 'antenna' that gives directivity.

you may also want to decompose the sum terms either

exp(j*x)=cos(x)+j*sin(x)

or

cos(x)=.5*(exp(j*x)+exp(-j*x))

sin(x)=-.5*j*(exp(j*x)-exp(-j*x))

but again, while the term exp(t) may try to pull up every thing, the tern exp(-n^2) is faster bringing it all down to zero.

[Walter Roberson]

You cannot calculate subtle variations in the values unless you use at least 192881416 digits in your calculation.

If you were to remove everything except the sum (and its multipliers) then you could take the log of the sum and get something marginally more tractable. It would still work out as just the first term but you would at least be able to represent it.

Now if you were to be varying t instead of r then... well, you are still pretty hopeless as long as D*t/a^2 is more than about 4.

But wait.. how do you know that the sin() should be over n*Pi*r/a and not extend to the exp() term, or not stop at the n ? You do not have any brackets that show the extent of the sin, so by convention it should extend as far as possible (or at the very least, to all terms multiplied together)