where is the error in my code?

Hi every one,

when I run the following code, I get a plot with unexpected two lines as circled by the red lines in the image attached! where is the error because we do not expect these lines to appear according to the physical equations?

V1 = @(r,w) -acosh(10*(w/(1600*r + 21))^(1/2))/20000000000
V2 = @(r,w) acosh(10*(w/(1600*r + 21))^(1/2))/20000000000
% Define function to be integrated
fun = @(x,r,w)0.0018./((w./((cosh(10^10.*x./0.5)).^2)-(r.*16+0.21)).^(1/2));
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).*log10((exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)))))./(r.*(integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)))))-(exp(-37.45.*r).*(70.31));
fimplicit(www,[0 5 0 0.075],'MeshDensity',500, 'LineWidth',1.5),grid

[Torsten]

I suggest you make a surface plot of www in the region of interest to see what happens.

help surf

[Abdallah Qaswal]

May you please help me plot the code with surf? because I tried using surf but I got an error, see here please:

V1 = @(r,w) -acosh(10*(w./(1600*r + 21))^(1/2))/20000000000
V2 = @(r,w) acosh(10*(w./(1600*r + 21))^(1/2))/20000000000
% Define function to be integrated
fun = @(x,r,w)0.0018./((w./((cosh(10^10.*x./0.5)).^2)-(r.*16+0.21)).^(1/2));
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))))/(r.*(integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)))))-(exp(-37.45.*r).*(70.31));
surf(www,[0 5 0 0.075],'MeshDensity',500, 'LineWidth',1.5),grid 

[Abdallah Qaswal]

when I delete the log part (log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)))), the straight linear line on the left is removed, so I think the problem is in writing the common logarithm, right?

Also, the horizontal line was removed when I change the range of R values from [0 0.075] to [0.01 0.075] so I removed the 0 to avoid infinity.

But the line that is remained is the left linear line! any idea to remove it?

[Torsten]

The part of your formula

...^(1/2))).log10(...

is incorrect. There must be a * or / before the log10.

w = 0:0.05:5;
r = 0.001:0.001:0.075;
[W,R]=meshgrid(w,r);
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))))/(r.*(integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)))))-(exp(-37.45.*r).*(70.31));
surf(W,R,www(W,R))

Dot indexing is not supported for variables of this type.

Error in solution (line 4)
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))))/(r.*(integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)))))-(exp(-37.45.*r).*(70.31));

[Abdallah Qaswal]

I edited the code according to your suggestions but unfortunately, the code still in error:

w = 0:0.05:5;
r = 0.001:0.001:0.075;
[W,R]=meshgrid(w,r);
V1 = @(r,w) -acosh(10*(w./(1600*r + 21)).^(1/2))/20000000000;
V2 = @(r,w) acosh(10*(w./(1600*r + 21)).^(1/2))/20000000000;
% Define function to be integrated
fun = @(x,r,w)0.0018./((w./((cosh(10^10.*x./0.5)).^2)-(r.*16+0.21)).^(1/2));
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).*log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))))./(r.*(arrayfun(@(r,w)integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)),R,W))))-(exp(-37.45.*r).*(70.31));
surf(W,R,www(W,R))

Warning: Inf or NaN value encountered.

[Abdallah Qaswal]

where is the problem?

[Torsten]

See the modified code from above.

You will have to check where and why the Inf and NaN values in the surface plotting appear.

I suggest you make a matrix M and inspect it:

w = 0:0.05:5;
r = 0.001:0.001:0.075;
[W,R]=meshgrid(w,r);
V1 = @(r,w) -acosh(10*(w./(1600*r + 21)).^(1/2))/20000000000;
V2 = @(r,w) acosh(10*(w./(1600*r + 21)).^(1/2))/20000000000;
% Define function to be integrated
fun = @(x,r,w)0.0018./((w./((cosh(10^10.*x./0.5)).^2)-(r.*16+0.21)).^(1/2));
www = @(w,r)5.124-6.4*10^-6.*(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))).*log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2))))./(r.*(arrayfun(@(r,w)integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)),R,W))))-(exp(-37.45.*r).*(70.31));
WWW = www(W,R);

Warning: Inf or NaN value encountered.

M(:,:,1) = W;
M(:,:,2) = R;
M(:,:,3) = WWW;
M(:,:,3)

ans = 75×101

1.0e+65 * -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000

[Abdallah Qaswal]

Thnak you! but this did not solve the original question because I would like to plot the relationship between R and W only without www, which is not a physical parameter in my case here!

As I told you, the issue is with the log part in the equation, any idea?

[Abdallah Qaswal]

Honestly, I do not understand the approach you suggested but what I noticed is that substituting the values on the undesired line on the left yields undefined values in the common logarithm like 0 or negative value, is there any function to exclude such undefined values?

[Abdallah Qaswal]

like just to take the real numbers

[Torsten]

Honestly, I do not understand the approach you suggested

MATLAB plots nonsense if www returns Inf or NaN.

So you will have to detect where this nonsense stems from in your definition of www.

[Abdallah Qaswal]

yes, the nonsens is (log of an exponential value), this one : log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)))) because once I remove it ( which does not affect the plot significantly), the plot was a desired one , do you have any idea how to address such log of an exponential function?

[Torsten]

log10(exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)))) =

-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)) / log(10)

Maybe this helps.

[Abdallah Qaswal]

Unfortunately, this also did NOT work, any other alternatives?

[Abdallah Qaswal]

I even put log 10 = 2.3, but still no desired results, here is the last modified code :

V1 = @(r,w) -acosh(10*(w/(1600*r + 21))^(1/2))/20000000000
V2 = @(r,w) acosh(10*(w/(1600*r + 21))^(1/2))/20000000000
% Define function to be integrated
fun = @(x,r,w)2.3*r.*0.0018./((w./((cosh(10^10.*x./0.5)).^2)-(r.*16+0.21)).^(1/2));
www = @(w,r)5.124+6.4*10^-6.*215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)).*exp(-215.6*0.5*(w.^(1/2)-(r.*16+0.21).^(1/2)))./(integral(@(x)fun(x,r,w),V1(r,w),V2(r,w)))-(exp(-37.45.*r).*(70.31));
fimplicit(www,[0 5 0.001 0.075],'MeshDensity',500, 'LineWidth',1.5),grid