If condition inside integration

Hi,

I am trying to keep if condition inside integration of exponential function and solve the integral in Matlab. c(x) is 4 at x =1, 4, 7, 10, 13 otherwise zero. Any help would be appreciated.

Answer 1

[Walter Roberson]

I am pretty sure you do not want that definition of c(x)

syms x t
c(x) = piecewise(x == 1 | x == 4 | x == 7 | x == 10 | x == 13, 4, 0)

c(x) = 

inner = c(x) * exp(-2*(t-x))

inner = 

y(t) = simplify(int(inner, x, 0, t))

y(t) = 

string(y)

ans = "piecewise(in(t, 'real'), 0)"

fplot(y, [-20 20])

This happens because your c(x) definition is discontinuous, and the width of the event x == 4 (or each of the other values) is 0, so the integral at those points is 0.

Compare:

c2(x) = (dirac(x-1) + dirac(x-4) + dirac(x-7) + dirac(x-10) + dirac(x-13)) * 4

c2(x) = 

inner2 = c2(x) * exp(-2*(t-x));
y2(t) = simplify(int(inner2, x, 0, t))

y2(t) = 

string(y2)

ans = "4*heaviside(t - 1)*exp(-2*t)*exp(2) + 4*heaviside(t - 4)*exp(-2*t)*exp(8) + 4*heaviside(t - 7)*exp(-2*t)*exp(14) + 4*heaviside(t - 10)*exp(-2*t)*exp(20) + 4*heaviside(t - 13)*exp(-2*t)*exp(26)"

fplot(y2, [-20 20])

This defines c(x) in terms of a distribution rather than in terms of points, giving meaning to the integral at those values.

[Dharma Khatiwada]

Thank you so much Walter! Your lower answer is something I was trying to do.

I need to insert the value of y2(t) in the loop below. How do I do that? lambda, b, d, e are all constant. Your suggestion will be greatly appreciated!

t=zeros(N,1);
V=zeros(N,1);
K=zeros(N,1);
%Run loop
% for j=1:N
t(1)=0;
V(1)=200; 
K(1)=625; 
for i=1:N
    for j=1:N
        t(i+1)=t(i)+dt;
        %Using Euler forward methodfor dV/dt
        V(i+1)=V(i)-lambda*dt*V(i)*log(V(i)/K(i)); 
        %Using Euler forward method for dK/dt
        K(i+1)=K(i)+b*dt*V(i)-d*dt*K(i)*(V(i).^(2/3))-e*dt*y2(t)*K(i);
    end
end

[Walter Roberson]

syms x t
c2(x) = (dirac(x-1) + dirac(x-4) + dirac(x-7) + dirac(x-10) + dirac(x-13)) * 4
inner2 = c2(x) * exp(-2*(t-x));
y2 = matlabFunction(simplify(int(inner2, x, 0, t)));


t=zeros(N,1);
V=zeros(N,1);
K=zeros(N,1);
%Run loop
% for j=1:N
t(1)=0;
V(1)=200; 
K(1)=625; 
for i=1:N
    for j=1:N
        t(i+1)=t(i)+dt;
        %Using Euler forward methodfor dV/dt
        V(i+1)=V(i)-lambda*dt*V(i)*log(V(i)/K(i)); 
        %Using Euler forward method for dK/dt
        K(i+1)=K(i)+b*dt*V(i)-d*dt*K(i)*(V(i).^(2/3))-e*dt*y2(t)*K(i);
    end
end

[Dharma Khatiwada]

Thank you again, Walter. Just a confusion with y2(t) in this line. It says 'left and right sides have different number of elements. Should I try y2(t(i))?

K(i+1)=K(i)+b*dt*V(i)-d*dt*K(i)*(V(i).^(2/3))-e*dt*y2(t)*K(i);

[Walter Roberson]

Yes, t(i) would make sense in that context.

[Dharma Khatiwada]

Hi Walter,

Just one more quick question with the value of c2(x) in the line:

c2(x) = (dirac(x-1) + dirac(x-4) + dirac(x-7) + dirac(x-10) + dirac(x-13)) * 4

Normally, we should expect dirac function giving the output i.e. infinity if for example, when x=1 otherwise zero. In my context, I am expecting that to be 4 i.e. c2(x)=4. I am not sure whether it is returning to that value. My concern is only on the value of c2(x).

I found something similar in other Matlab posting. Please suggest if I also need to normalize like below.

Thank you again!

x = -1:0.1:1;
y = dirac(x);
idx = y == Inf; % find Inf
y(idx) = 1;     % set Inf to finite value

[Walter Roberson]

dirac(0)

ans = Inf

syms x
int(dirac(x), x, -1, 1)

ans = 

1

The Dirac Delta is not really a function. dirac(0) is not really inf . dirac() is defined such that the integral across 0 is 1. So what happens with int() of dirac is correct, and the Inf is not really correct.

There are different ways of defining Dirac. One of the ways is as the limit of a rectangle n units high and 1/n wide, as n approaches infinity: the area is fixed, but as the width approaches 0, the height approaches infinity. Saying that it is infinity such as dirac(0) shows, is a short-hand that is not really true.