Evolution of a parameter over time

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

Answer 1

[Ameer Hamza]

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.