why power series plot is not matching with ode

1 view (last 30 days)
t=1:20;
a(1)=1;
b(1)=1;
c(1)=1;
c1=10;
c2=28;
c3=8/3;
sumx=0;
sumy=0;
sumz=0;
for i=1:10
a(i+1)=(c1*(b(i)-a(i)))/(i+1);
b(i+1)=((c2-c(i))*a(i)-b(i))/(i+1);
c(i+1)=(a(i)*b(i)-c3*c(i))/(i+1);
sumx=sumx+a(i)*t^i;
sumy=sumy+b(i)*t^i;
sumz=sumy+c(i)*t^i;
x=sumx;
y=sumy;
z=sumz;
end
plot(t,x) not matching
  8 Comments
shiv gaur
shiv gaur on 22 Feb 2022
Edited: Walter Roberson on 22 Feb 2022
t=1:20;
a(1)=1;
b(1)=1;
c(1)=1;
c1=10;
c2=28;
c3=8/3;
sumx=0;
sumy=0;
sumz=0;
for i=1:10
a(i+1)=(c1*(b(i)-a(i)))/(i+1);
b(i+1)=((c2-c(i))*a(i)-b(i))/(i+1);
c(i+1)=(a(i)*b(i)-c3*c(i))/(i+1);
sumx=sumx+a(i)*t^i;
sumy=sumy+b(i)*t^i;
sumz=sumy+c(i)*t^i;
x=sumx;
y=sumy;
z=sumz;
end
program of power series pl help to correct graph
Torsten
Torsten on 22 Feb 2022
According to
the coefficients of the power series for x(t), y(t) and z(t) of the Lorenz system are much more difficult than
a(i+1)=(c1*(b(i)-a(i)))/(i+1);
b(i+1)=((c2-c(i))*a(i)-b(i))/(i+1);
c(i+1)=(a(i)*b(i)-c3*c(i))/(i+1);

Sign in to comment.

Answers (2)

Walter Roberson
Walter Roberson on 22 Feb 2022
syms x
f = log(x+1)
f = 
xtay = taylor(f, x, 0, 'order', 10)
xtay = 
subplot(2,1,1); fplot(f, [-30 30])
subplot(2,1,2); fplot(xtay, [-30 30])
Just because you take a truncated taylor series doesn't mean that the truncated taylor series is at all accurate as you get away from the expansion point.
  4 Comments
Walter Roberson
Walter Roberson on 23 Feb 2022
You have the complicated plot from https://www.mathworks.com/matlabcentral/answers/1655855-why-power-series-plot-is-not-matching-with-ode#comment_2000660 . If you count, it has about 50 inflection points -- places where the plot changes direction between up and down. That means that it has at least 50 places where the derivative is 0. But in order to have 50 places with derivative 0 expressed as a polynomial (which your power series expansion is), you need a polynomial of degree 50; a polynomial of degree 10 has no chance.
Walter Roberson
Walter Roberson on 23 Feb 2022
Also, you are plotting only 20 points; you cannot possibly see 50 changes of direction in only 20 points. You would need at least 101 points to see 50 changes of direction.
If you try to extend numerically to higher degrees, higher maximum i, then i = 15 is as high as your code can go before some of the coefficients overflow to infinity.
If you switch to symbolic computation and ask to go higher degree.. then it takes a long time. And you tend to get points evaluated to be on the order of 10^27000.
At this point you have several possibilities to consider:
  • that if you evaluated to a "high enough" degree symbolically, that you just might get a decent interpolation over t = 1 to t = 20
  • that, because the minimum degree for replicating the plot is at least 50, and you are evaluating out to t = 20, you are going to have constant * 20^50 ~~ 10^65 and you are hoping for final results in the order of +/- 20... you are going to need to evaluate symbolically to extra digits or else the values you want will disappear as noise
  • that since the system appears to have an infinite number of changes of direction (we can mentally interpolate that there will be many more peaks at increasing x), that any finite polynomial approximation might perhaps be too badly conditioned to be useful, calculations that might work in theory but are more like balancing on the edge of a pin numerically
  • your implementation for calculating the coefficients might be wrong. But as I have taken the equations and implemented them in a different software package, I know that even if you got a detail wrong, that the proper plot looks a lot like the first plot, and an order 10 series expansion for x looks a lot like what you plotted, so even if you got a detail wrong, the exact same challenges would be sure to happen even if you had the details perfect.

Sign in to comment.


shiv gaur
shiv gaur on 23 Feb 2022
dear sir this is the paper which is solve by power series method
  12 Comments
shiv gaur
shiv gaur on 24 Feb 2022
pl solve the lorenz system by power series
Walter Roberson
Walter Roberson on 24 Feb 2022
I have been working for about 18 hours. I need rest.

Sign in to comment.

Categories

Find more on Mathematics in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!