0

I am calculating the solution of Chen's chaotic system using differential transform method. The code that I am using is:

x=zeros(1,7);
x(1)=-0.1;
y=zeros(1,7);
y(1)=0.5;
z=zeros(1,7);
z(1)=-0.6;
for k=0:5
    x(k+2)=(40*gamma(1+k)/gamma(2+k))*(y(k+1)-x(k+1));
    sum=0;
    for l=0:k
        sum=sum+x(l+1)*z(k+1-l);
    end
    y(k+2)=(gamma(1+k)/gamma(2+k))*(-12*x(k+1)-sum+28*y(k+1));
    sum=0;
    for l=0:k
        sum=sum+x(l+1)*y(k+1-l);
    end
    z(k+2)=(gamma(1+k)/(1+k))*(sum-3*z(k+1));
end
s=fliplr(x);
t=0:0.05:2;
a=polyval(s,t);
plot(t,a)

What this code does is calculate x(k), y(k) and z(k) these are the coefficients of the polynomial that is approximating the solution. The solution x(t) = sum_0^infinity x(k)t^k, and similarly the others. But this code doesn't give the desired output of a chaotic sequence the graph of x(t) that I am getting is:

enter image description here

EBH
  • 10,350
  • 3
  • 34
  • 59
user311790
  • 192
  • 10

1 Answers1

1

This is not an answer, but a clearer and more correct (programmatically speaking) to write your loop:

for k = 1:6
    x(k+1)=(40*1/k)*(y(k)-x(k));
    temp_sum = sum(x(1:k).*z(k:-1:1),2);
    y(k+1) = (1/k)*(-12*x(k)-temp_sum+28*y(k));
    temp_sum = sum(x(1:k).*y(k:-1:1),2);
    z(k+1) = (1/k)*(temp_sum-3*z(k));
end

The most important issue here is not overloading the built-in function sum (I replaced it with temp_sum. Other things include vectorization of the inner loops (using sum...), indexing that starts in 1 (instead of writing k+1 all the time), and removing unnecessary calls to gamma (gamma(k)/gamma(k+1) = 1/k).

EBH
  • 10,350
  • 3
  • 34
  • 59
  • wouldn't there be an array indexing error in line 3 of your code when k=1; and also gamma(1+k)/gamma(2+k)=1/(k+1) – user311790 May 24 '17 at 17:23
  • @Upstart, it doesn't throw an error, but it was wrong, I corrected it now. And also what I meant to write was `gamma(k)/gamma(k+1) = 1/k`. – EBH May 24 '17 at 17:40
  • it still throws an error in matlab.can you explain the third line to me? – user311790 May 24 '17 at 17:44
  • I made it a little clearer. It says: take the elements from 1 to k in `x`, and multiply element-wise (`.*`) by elements k to 1 in `z`. Then sum all the results and assign to `temp_sum`. BTW, I have checked it and it gives the exact same result with no errors, make sure you type `clear` before you run this. – EBH May 24 '17 at 17:50
  • what does that ",2" mean? and also when u calculate z why have you used gamma(k)/k? – user311790 May 24 '17 at 17:52
  • It's optional, it means "sum over the second dimension" (i.e the columns). If you omit it and the input is a vector, then Matlab always sum it along the non-singleton dimension. – EBH May 24 '17 at 17:56
  • but still this is not the answer to my question – user311790 May 24 '17 at 17:58
  • You have to remember the `k` goes now from 1 to 6 (and not 0 to 5), so instead of taking `gamma(1+k)/(1+k)` we need to take `gamma(k)/k` – EBH May 24 '17 at 17:58
  • I know it's not an answer, that the first thing I wrote in this non-answer... but writing a clear and efficient code is also important, and reduce the probability of mistakes (and improve the readability of the code and the question). – EBH May 24 '17 at 17:59
  • in my code there was a typo in the z calculation also there was gamma(1+k)/gamma(2+k) – user311790 May 24 '17 at 18:03
  • how do i write $\sum_{k=0}^10 X(k)(t-t_0)^{ka}$ in vectorize? – user311790 Jul 31 '17 at 10:47
  • @Upstart please ask new questions in new posts, not in comments. Also, don't use latex here, it isn't interpreted. – EBH Jul 31 '17 at 10:49