MATLAB: Evolution of a parameter over time

codedifferential equationserrorgraphiterationMATLABtime series

I am elaborating a code to simulate the evolution of a specific parameter (D) over time.
This my code:
b= 70;
r= 2.78;
u = 0.4;
da = 0.12;
sigma1 = 26; % MPa


sigma2 = 1; % MPa
sigma3 = 1; % MPa
t(1)=0
D(1) = 1e-5;
%Calculation
for i=1:400
t(i+1)= t(i)+0.01
H(i)= D(i)-da;
%Switch Function
if H(i)>=0
    s(i)=H(i);
else
    s(i)=0;
end
sigma = sigma1-sigma3;
sigmam = (sigma1+sigma2+sigma3)/3;
sigmadano = sigma*((2/3)*(1+u)+3*(1-2*u)*((sigmam/sigma)^2))^0.5
D(i+1) = ((sigmadano/(b*(1-D(i))*(1-s(i))))^r)*t(i+1));
%Set D maximum
if D(i+1)>1
    D(i+1)=1
break
end
end
The first value of D is 1e-5 and it was considered as D(1)=1e-5 in my code. Time is "t" in my code and t(1)=0
My doubt is how can I represent the evolution of D parameter using MATLAB. I think it would be something like D(i+1), as it is shown in my code, but something is wrong.
I got this graph:
But the correct is:
I don't know where is my error.
Thanks

Best Answer

Since you are trying to numerically approximate a derivative, you also need to add all the previous values weighted with time. Replace the following line in your code
D(i+1) = ((sigmadano/(b*(1-D(i))*(1-s(i))))^r)*t(i+1));
with
D(i+1) = D(i) + ((sigmadano/(b*(1-D(i))*(1-s(i))))^r)*(t(i+1)-t(i));
Here I am using left point rule which is quite simple but is giving good estimation in your case.
Related Question