I wrote the code but still the answer is not correct using delta function

the two recursive formula are as follows:

C(k,h+1)=D1*(k+2)*(k+1)*C1(k+2,h)+lambda*omega*dirac delta(k,h)-lambda C1(k,h)-k2*((Sum(r=0 to k) sum(s-0 to h)C1(k-r,s)*C2(r,h-s)

2nd equation

C2(k,h+1)=D2*(k+2)*(k+1)C2(k+2,h)-K2*((sum(r=0 to k) sum(s=0 to h) C1(k-r,h)*C2(r,h-s)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D_1=0.0001
D_2=0.0001
lambda=3;
omega=1
k_2=1

syms x
syms t

C1=zeros(200,200);
C2=zeros(200,200);

series(x,t)=sym(zeros(1,1));

for k=2:1:150
   C1(1,1)=1;
   C1(k,1)=0;

end
for k=2:1:18
  C2(1,1)=6;  
  C2(k,1)=0;
end

for h=1:1:150
    for k=1:1:10
        S=0;
        for r=1:1:k
            for s=1:1:h
             if r==1&& h+1==1 for delta condition
            D=1;
        else
            D=0;
        end    

              S=S+C1(k-r+1,s)*C2(r,h-s+1);    
              end
              end

          C1(k,h+1)=(D_1*k*(k+1)*C1(k+2,h))+(lambda*omega*D)-(lambda*C1(k,h))-(k_2*S)/h;
          C2(k,h+1)=((D_2*k*(k+1)*C2(k+2,h))-(k_2*S))/h;
      end
  end

Answer 1

[Walter Roberson]

Your formula has dirac delta(k,h) which would typically indicate the k'th derivative of the dirac delta with respect to h.

The 0'th derivative or normal dirac delta, dirac(h), is a distribution (not a function) that is nominally defined as +infinity at h == 0 and as 0 otherwise, but is used more for its property that constant * dirac(h) integrates to the given constant at the point where h is 0 (this is something that cannot happen with a true mathematical function; it takes a distribution.) The k'th derivative of the dirac delta is typically (-1)^k * infinity when evaluated in isolation.

Basically a dirac delta only make sense to use if you are going to do symbolic integration (not numeric integration either): a dirac delta that is not being integrated is going to contribute either 0 or an infinity to the function and that infinity is going to mess up your calculations.

Your code has

        if r==1&& h+1==1 for delta condition
            D=1;
        else
            D=0;
        end

which attempts to define D as 0 except where r == 1 and h == 0, and r is (for reasons that are not immediately clear) looping to k. That is not correct outside of an integral: outside of an integral, D should become +/- infinity (the sign depending on r).

[YOGESHWARI PATEL]

r depends on k.so depending upon value of k the loop will run for summation.If this is not correct then can u suggest that with term i should replace it so my answer is correct

[Walter Roberson]

Change

if r==1&& h+1==1 for delta condition
          D=1;
      else
          D=0;
      end

to

if r==1&& h+1==1 for delta condition
          D=-inf;
      else
          D=0;
      end

Better yet, rewrite your summations so that each of the possibilities is summed separately. That would include not using D at all, replacing

(lambda*omega*D)

with

(-1)^.r * inf

in the cases where h is 0, and dropping the (lambda*omega*D) term completely for cases where h is non-zero.