Why my graph not same as research paper?

Question

nur 2022-6-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1736375-why-my-graph-not-same-as-research-paper

评论： Asif Solanki 2022-8-25

采纳的回答： Torsten

在 MATLAB Online 中打开

Hi, i want to ask why i did not get same graph as above?

This is my code that i have try:

%for t1=5%

x=0:0.2:10;

m=0.1;

C_0=1.0;

u_0=0.25;

D_0=0.45;

w_0=0.001;

p_0=0.02;

t1=5;

a=((((u_0)^2)/(4*D_0))+p_0);

b=((u_0)^2)/(4*D_0);

d=(1/(sqrt(2)*m));

e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t1/2))));

e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t1/2)));

T=d*(e1-e2);

f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));

k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));

l=(((u_0)*x)/(D_0));

n=exp(h);

o=erfc(f);

q=exp(i);

r=erfc(g);

s=exp(-p_0*T);

v=erfc(j);

y=erfc(k);

z=exp(l);

F=(((1/2).*n.*o)+((1/2).*q.*r));

G=(s.*(1-(1/2).*v-(1/2).*z.*y));

C1=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);

plot(x,C1)

hold on

%for t2=10%

x=0:0.2:10;

m=0.1;

C_0=1.0;

u_0=0.25;

D_0=0.45;

w_0=0.001;

p_0=0.02;

t2=10;

a=(((u_0)^2)/(4*D_0)+p_0);

b=((u_0)^2)/(4*D_0);

d=(1/(sqrt(2)*m));

e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t2/2))));

e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t2/2)));

T=d*(e1-e2);

f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));

k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));

l=(((u_0)*x)/(D_0));

n=exp(h);

o=erfc(f);

q=exp(i);

r=erfc(g);

s=exp(-p_0*T);

v=erfc(j);

y=erfc(k);

z=exp(l);

F=(((1/2).*n.*o)+((1/2).*q.*r));

G=(s.*(1-(1/2).*v-(1/2).*z.*y));

C2=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);

plot(x,C2)

hold off

hold on

%for t3=15%

x=0:10;

m=0.1;

C_0=1.0;

u_0=0.25;

D_0=0.45;

w_0=0.001;

p_0=0.02;

t3=15;

a=(((u_0)^2)/(4*D_0)+p_0);

b=((u_0)^2)/(4*D_0);

d=(1/(sqrt(2)*m));

e1=(log(2+sqrt(2)+(4+3*(sqrt(2)))*(tanh(m*t3/2))));

e2=(log(2+sqrt(2)-sqrt(2)*tanh(m*t3/2)));

T=d*(e1-e2);

f=((x-((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

g=((x+((((u_0)^2)+4*p_0*D_0)^1/2)*T)/(2*sqrt((D_0)*T)));

h=((((u_0)-((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

i=((((u_0)+((((u_0)^2)+4*p_0*D_0)^1/2))*x)/(2*D_0));

j=((x-(u_0)*T)/(2*sqrt((D_0)*T)));

k=((x+(u_0)*T)/(2*sqrt((D_0)*T)));

l=(((u_0)*x)/(D_0));

n=exp(h);

o=erfc(f);

q=exp(i);

r=erfc(g);

s=exp(-p_0*T);

v=erfc(j);

y=erfc(k);

z=exp(l);

F=(((1/2).*n.*o)+((1/2).*q.*r));

G=(s.*(1-(1/2).*v-(1/2).*z.*y));

C3=(((w_0)/(p_0))+(C_0-((w_0)/(p_0))))*F-(((w_0)/(p_0))*G);

plot(x,C3)

hold on

legend({'t=5','t=10','t=15'},'Location','southwest')

xlabel('Distance,x (meter)')

ylabel('Concentration,C(x,t)')

%ylim([0, 1])

xlim([0, 10])

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2022-6-8

The only thing that helps:

Compare the formulas from the article and yours.

Torsten 2022-6-8

@nur

I corrected the three errors in your code detected by @SALAH ALRABEEI

Looks better now, but concentrations still become negative. So your search must go on.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-6-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1736375-why-my-graph-not-same-as-research-paper#answer_981715

编辑：Torsten 2022-6-8

在 MATLAB Online 中打开

Sometimes it's better to start anew:

x=(0:0.1:10).';

C0 = 1.0;

tt = [5.0,10.0,15,0];

u0 = 0.25;

D0 = 0.45;

m = 0.1;

gamma0 = 0.02;

mu0 = 0.001;

C = zeros(numel(x),numel(tt));

for i=1:numel(tt)

t = tt(i);

T = 1/(sqrt(2)*m)*(log(2+sqrt(2)+(4+3*sqrt(2))*tanh(m*t/2))-...

log(2+sqrt(2)-sqrt(2)*tanh(m*t/2)));

F = 0.5*exp((u0-sqrt(u0^2+4*gamma0*D0))*x/(2*D0)).*...

erfc((x-sqrt(u0^2+4*gamma0*D0)*T)./(2*sqrt(D0*T)))+...

0.5*exp((u0+sqrt(u0^2+4*gamma0*D0))*x/(2*D0)).*...

erfc((x+sqrt(u0^2+4*gamma0*D0)*T)./(2*sqrt(D0*T)));

G = exp(-gamma0*T)*(1-0.5*erfc((x-u0*T)./(2*sqrt(D0*T)))-...

0.5*exp(u0*x/D0).*erfc((x+u0*T)./(2*sqrt(D0*T))));

C(:,i) = mu0/gamma0 + (C0-mu0/gamma0)*F - mu0/gamma0*G;

end

plot(x,C)

2 个评论
显示无隐藏无

nur 2022-6-9

thank you so much!!! this graph exactly same as in the research paper. :)

Asif Solanki 2022-8-25

Dear Sir,

My name is Asif Solanki I am a student persuing master's program in Italy. I am currently working on a small project where I do need to find the concentration of the pollutants using the Bear analytical method.

Therefore, would you like to assist me to find the solution?

Thank you very much

Best regards

Asif Solanki

请先登录，再进行评论。

Answer 2

Sam Chak 2022-6-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1736375-why-my-graph-not-same-as-research-paper#answer_981515

Hi @nur

One of the ways to find out is to determine the equilibrium points from the advection-dispersion equation.

I haven't checked the long equations. Can you verify if you plotted C or

? Sometimes, the transformations can be a little tricky.

2 个评论
显示无隐藏无

nur 2022-6-8

I plotted C. I hope you can help me ;(

Sam Chak 2022-6-8

编辑：Sam Chak 2022-6-8

Hi @nur

Try to make a little bit of value-added effort to the problem.

Another way is to plot the concentration C from the numerical solution of the advection-dispersion differential equation.

This way you can compare with the analytical solution obtained from the paper. Bear in mind that sometimes misprints can occur due to authors' mistake, or the production crew's mistake. So, what was shown on the paper might not be truly the analytical solution. Therefore, yes you have to check.

But I'd suggest you to take numerical solution approach because that is directly from the governing advection-dispersion law. The analytical solution, which is "human-processed equation" from the law.

Edit: Hey, check this out:

https://nnesmo2019.webspace.durham.ac.uk/practicals/session42/

One of the objectives is to Solve the one-dimensional advection dispersion equation using ode45 or ode15s.

请先登录，再进行评论。