How to write integral equation ?

I like to check if my code represent the equation. Please help

[David Hill]

What is your function D(c)? Or is it an array? Depending on what D(c) is, and assuming you want a numerical solution, it should not be a hard problem.

[Michael Crown]

D(c)

[Michael Crown]

It should return a single value. Not an array

[David Hill]

What is your function D with respect to c? For example:

D=@(c)5*c.^2+4*c-10;

[Michael Crown]

That is why I need help. If I have your email address, I can send the entire code to you. I am a novice with this. Thanks.

[Michael Crown]

My email: [e-mail] You can send me a mail and I will reply with the code if you aren’t comfortable sharing email here

[Walter Roberson]

It looks to me as if you should be using trapz()

[Michael Crown]

If you don’t mind, can I send the full code? I need help.

[David Hill]

Yes, attach full code and ask a specific question.

[Michael Crown]

attached @David Hill

All i need is to calculate DaveA and DaveB using the attached equation. I am not sure what i have i right.

Answer 1

[David Hill]

Not sure if your Ca is sorted. Yes, you have a stream of data that is being integrated. The proper function is trapz(). Providing your full code helped significantly. I assume you want your data sorted before integrating.

[Ca,idx]=sort(dataB(:,1));
Da=dataB(:,2);
Da=Da(idx);%aligns Da to sorted Ca
DaveA=trapz(Ca,Da)/(Ca(end)-Ca(1));

[David Hill]

The above is all you need if you have dataB. It will calculate DaveA per the equation you provided.

[Michael Crown]

Thanks so much!

[David Hill]

I do not know what matanoIndex is, so you will have to code that.

time = 501*60*60;
[pos,idx]=sort(data(:,1).*3e2);
c=data(:,2);
c=d(idx);
c=smooth(d,.2);
cMax=c(end);
cMin=c(1);
pos=smooth(c,.2);
matanoPos= pos(matanoIndex);
posn =(pos - matanoPos);
slope = 1/2*diff(posn(1:end-1))./diff(c(1:end-1))+diff(posn(2:end))./diff(c(2:end));	
i= cumsum(1/2*movsum(posn(2:end),2,'Endpoints','discard').*diff(c(1:end-1)));
c=c(2:end-1);
d=-1/(2*time).*slope.*i;
DaveM=trapz(c,d)/(cMax-cMin);

[Michael Crown]

Please check, D was calculated in the code and not given in input file like DA. I get error code: “x must be a vector” when I used this code for DaveM

[Michael Crown]

I realised that it retuned an answer without error now but it messed up my plot for D vrs C

[David Hill]

You will have to figure it out from here. conc was not needed (only the transpose of c). d is calculated in your code. I vectorized and eliminated your loops. Keep in mind, you smoothed for DaveM but did not smooth for DaveA. Also cMax and cMin were calculated before the elimination of the two end-points for the integration. I would suggest determining cMax and cMin just before integration.

[Michael Crown]

Is there another way in such that DaveM is calculated from the value I got for D vz c? This code you shared messed the plots.

[Michael Crown]

Hello, if you don’t mind, help modify code such that DaveM will be calculated from the already obtained D vrs c from my data without sorting . The one you sent to me isn’t working.

Answer 2

[Michael Crown]

Forgive me for too many requests. Please can you put the modification in the code and send back. I am not clear

[Michael Crown]

Hello David, I am not able to calcualte DaveM. can you help check?

Answer 3

[Michael Crown]

Now, i am totally confused, it didnt run

time = 223*60*60;	% diffusion time in seconds	
data = load ('File 223.txt'); % read in data
dataB = load ('File 2233.txt');

[pos,idx]=sort(data(:,1)); 

c=data(:,2 ); 
D= D(idx ); 
c=smooth(c,.2 ); 
D = smooth(D, .2);
cMax=c(end ); 
cMin=c(1 ); 
pos=smooth(c,.2 ); 


matanoPos= pos(matanoIndex ); 
posn =(pos - matanoPos ); 
slope = 1/2*diff(posn(1:end-1))./diff(c(1:end-1))+diff(posn(2:end))./diff(c(2:end )); 
i= cumsum(1/2*movsum(posn(2:end),2,'Endpoints','discard').*diff(c(1:end-1 ))); 
c=c(2:end-1 ); 
D=-1/(2*time).*slope.*i ; 
DaveM=trapz(c,d)/(cMax-cMin);


[cA,idx] = sort(dataB(:,1));
DA=DA(idx);
DaveA = trapz(cA,DA)/(cAMax-cAMin); %cAMax = cA(end), cAMin = cA(1)

semilogy (c,D,'r*') %create a plot using a base 10 logarithmic scale for y-axis and a linear scale for x-axis
hold on
semilogy(cA,DA, 'k*')


title('D of the Cu–Zn binary system at 700degC, 501hrs')
xlabel('Concentration, Zn(at.%)')
ylabel('Interdiffusion Coefficient (m2/s)')
legend('Mike','','Author')

Answer 4

[Michael Crown]

Hello Walter Roberson,

I still need your help on my codes. Attahced codes; mainBMFull and mainSF are written for two methods of calcualting 'D'

My mainSF code works well and gives me exactly what i want. However, mainBMFull is not accurate. Please, can you help check this?

let me know what else you need to help me.