taylor series method expansion

I wrote the following code

function yb=taylor(f,a,b,ya,n)
syms x y;
h=(b-a)/n;
y(1)=ya;
y(2)=f;
ht=h.^(0:5)./factorial(0:5)
for i=2:3
      y(i+1)=diff(y(i),x);
end
y

I got the output as follows

>> f=@(x)x^2+y^2;
>> taylor(f,0,0.8,1.5,4)
ht =
    1.0000    0.2000    0.0200    0.0013    0.0001    0.0000
y =
[ 3/2, x^2 + y(x)^2, 2*x + 2*y(x)*diff(y(x), x), 2*diff(y(x), x)^2 + 2*y(x)*diff(y(x), x, x) + 2]

Suggest me how to perform the following actions

1) I want to evaluate y at different values of x

2) Multiply ht and y elements and sum

[PJS KUMAR]

please suggest me the code in MATLAB to expand the taylor series method for the following problem

y' = x^2+y^2, for the values of x = 1 (0.5) 3 with y(1)=0.8

[Torsten]

function main
  x = linspace(1,3,5);
  y = zeros(numel(x))
  deltax = x(2)-x(1);
  y(1) = 0.8;
  for i = 1:numel(x)-1
    D = dy(x(i),y(i));
    y(i+1) = y(i)+deltax*D(1)+deltax^2/2*D(2)+deltax^3/6*D(3)+deltax^4/24*D(4)+deltax^5/120*D(5);
  end
plot(x,y)
end

function D = dy(x,y)
  f = x^2+y^2;
  df = 2*x+2*y*f;
  d2f = 2+2*(f^2+y*df);
  d3f = 2*(2*f*df+f*df+y*d2f);
  d4f = 2*(2*(df^2+f*d2f)+df^2+f*d2f+f*d2f+y*d3f);
  D = [f df d2f d3f d4f];
end

[PJS KUMAR]

I wrote the following code for Taylor series expansion

function yb=taylor(f,a,b,ya,n)
syms x y;
h=(b-a)/n;
y(1)=ya;
y(2)=f;
ht=h.^(0:3)./factorial(0:3)
for i=2:3
      y(i+1)=diff(y(i),x);
end
s=sum(ht.*y)

I run code and got the following output

>> syms x y(x);
>> f=@(x)x^2+y^2;
>> taylor(f,1,3,0.8,4)
s =
x/4 + (y(x)*diff(y(x), x))/4 + diff(y(x), x)^2/24 + y(x)^2/2 + (y(x)*diff(y(x), x, x))/24 + x^2/2 + 101/120

Now i want to evaluate the output 's' at different values of x and y, where x is from a to b and y(a)=0.8

[Walter Roberson]

y' = x^2+y^2, for the values of x = 1 (0.5) 3 with y(1)=0.8

I am having difficulty understanding the initial conditions. Could you confirm that you are referring to

syms y(x)
>> simplify(dsolve(diff(y,x) == x^2+y^2, y(1)==0.8))
ans =
piecewise(C5 ~= 0, (x*(4*besselj(-1/4, 1/2)*besselj(-3/4, x^2/2) + 4*besselj(1/4, 1/2)*besselj(3/4, x^2/2) + 5*besselj(-3/4, 1/2)*besselj(3/4, x^2/2) - 5*besselj(3/4, 1/2)*besselj(-3/4, x^2/2)))/(4*besselj(1/4, 1/2)*besselj(-1/4, x^2/2) - 4*besselj(-1/4, 1/2)*besselj(1/4, x^2/2) + 5*besselj(-3/4, 1/2)*besselj(-1/4, x^2/2) + 5*besselj(3/4, 1/2)*besselj(1/4, x^2/2)))

and you want taylor expansion of that without having solved the differential equation ?

Answer 1

[PJS KUMAR]

function D = dy(x,y)
  f = x^2+y^2;
  df = 2*x+2*y*f;
  d2f = 2+2*(f^2+y*df);
  d3f = 2*(2*f*df+f*df+y*d2f);
  d4f = 2*(2*(df^2+f*d2f)+df^2+f*d2f+f*d2f+y*d3f);
  D = [f df d2f d3f d4f];
end

suggest me to change the above code as given below, so that the differentiation is calculated by using diff function, without giving the differentiation directly in the program.

function D = dy(x,y)
  f = x^2+y^2;
  df = diff(f);
  d2f = diff(df);
  ---
  ---

[Torsten]

Just for testing: Does this work ?

If it works, you should do the symbolic differentiations once and pass the obtained symbolic expressions f, df, d2f, d3f, d4f to "dy" during the numerical integration.

function D = dy(xnum,ynum)
  syms x y(x)
  f = x^2+y^2;
  df = diff(f,x);
  d2f = diff(df,x);
  d3f = diff(d2f,x);
  d4f = diff(d3f,x);
  fnum = double(subs(f,[x,y],[xnum,ynum]));
  dfnum = double(subs(df,[x,y],[xnum,ynum]));
  d2fnum = double(subs(d2f,[x,y],[xnum,ynum]));
  d3fnum = double(subs(d3f,[x,y],[xnum,ynum]));
  d4fnum = double(subs(d4f,[x,y],[xnum,ynum]));
  D = [fnum,dfnum,d2fnum,d3fnum,d4fnum];
end

[PJS KUMAR]

I wrote the program as below

function main
 x = linspace(1,3,5);
 % y = zeros(numel(x));
 h = x(2)-x(1);
 y(1) = 0.8;
 for i = 1:numel(x)-1
   D = dy(x(i),y(i));
    y(i+1) = y(i)+h*D(1)+h^2/2*D(2)+h^3/6*D(3)+h^4/24*D(4)+h^5/120*D(5);
 end
[x' y'] 
end

function D = dy(xnum,ynum)
  syms x y(x)
  f = x^2+y^2;
  df = diff(f,x);
  d2f = diff(df,x);
  d3f = diff(d2f,x);
  d4f = diff(d3f,x);
  fnum = double(subs(f,[x,y],[xnum,ynum]));
  dfnum = double(subs(df,[x,y],[xnum,ynum]));
  d2fnum = double(subs(d2f,[x,y],[xnum,ynum]));
  d3fnum = double(subs(d3f,[x,y],[xnum,ynum]));
  d4fnum = double(subs(d4f,[x,y],[xnum,ynum]));
  D = [fnum,dfnum,d2fnum,d3fnum,d4fnum];
end

It is giving error as

>> main
Index exceeds matrix dimensions.

Error in sym/subs>normalize (line 193)
        Y{i} = fun2sym(Y{i},argnames(X{i}));

Error in sym/subs>mupadsubs (line 147)
[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 135)
    G = mupadsubs(F,X,Y);

Error in main>dy (line 20)
  fnum = double(subs(f,[x,y],[xnum,ynum]));

Error in main (line 7)
   D = dy(x(i),y(i));

kindly correct the program

[Torsten]

I don't have the Symbolic Toolbox - so I can't make a test.

Try

function D = dy(xnum,ynum)
syms x y(x)
f = x^2+y^2;
df = diff(f,x);
df = subs(df,diff(y(x),x),f);
d2f = diff(df,x);
d2f = subs(d2f,diff(y(x),x),f);
d3f = diff(d2f,x);
d3f = subs(d3f,diff(y(x),x),f);
d4f = diff(d3f,x);
d4f = subs(d4f,diff(y(x),x),f);
fnum = double(subs(f,[x,y],[xnum,ynum]));
dfnum = double(subs(df,[x,y],[xnum,ynum]));
d2fnum = double(subs(d2f,[x,y],[xnum,ynum]));
d3fnum = double(subs(d3f,[x,y],[xnum,ynum]));
d4fnum = double(subs(d4f,[x,y],[xnum,ynum]));
D = [fnum,dfnum,d2fnum,d3fnum,d4fnum];
end

[PJS KUMAR]

It is working correctly. A small modification needed in 'subs' syntax

fnum = double(subs(f,{x,y},{xnum,ynum}));
dfnum = double(subs(df,{x,y},{xnum,ynum}));
d2fnum = double(subs(d2f,{x,y},{xnum,ynum}));
d3fnum = double(subs(d3f,{x,y},{xnum,ynum}));
d4fnum = double(subs(d4f,{x,y},{xnum,ynum}));

[PJS KUMAR]

1) We used 'D' function for 4 derivatives. suggest the code to expand this D function for 'n' number of derivatives.

2) Modify the code so that the function shall be given from the command prompt

[Walter Roberson]

"suggest the code to expand this D function for 'n' number of derivatives."

Loop. Store the derivatives in a symbolic array.

You would not need to loop doing the subs(): you could subs() on the entire array at one time.

You can vectorize the calculation of the factorials,

factorial(0:N)

You can vectorize the calculation of the powers of h:

h.^(0:N)

[Torsten]

function D = dy(n,xnum,ynum)
 syms x y(x)
 dfstart= x^2+y^2;
 dfold = dfstart;
 D(1) = double(subs(dfold,{x,y},{xnum,ynum}));
 for i=2:n
   dfnew = subs(diff(dfold,x),diff(y(x),x),dfstart);
   D(i) =  double(subs(dfnew,{x,y},{xnum,ynum}));
   dfold = dfnew;
 end
end

I leave it to you to modify the main program appropriately.

Best wishes

Torsten.

[PJS KUMAR]

The modified 'D' function for order 'n' is working.

1) suggest me to write the single function by combining both 'main' and 'D' functions.

2) The function shall be called with the parameters f,xa,xb,ya,h

function main
 x = linspace(1,3,5);
 h = x(2)-x(1);
 y(1) = 0.8;
 for i = 1:numel(x)-1
   D = dy(7,x(i),y(i));
    y(i+1) = y(i)+h*D(1)+h^2/2*D(2)+h^3/6*D(3)+h^4/24*D(4)+h^5/120*D(5);
 end
[x' y'] 
end

function D = dy(n,xnum,ynum)
 syms x y(x)
 dfstart= x^2+y^2;
 dfold = dfstart;
 D(1) = double(subs(dfold,{x,y},{xnum,ynum}));
 for i=2:n
   dfnew = subs(diff(dfold,x),diff(y(x),x),dfstart);
   D(i) =  double(subs(dfnew,{x,y},{xnum,ynum}));
   dfold = dfnew;
 end
end

[Torsten]

Now the meal is prepared. You only have to put it on the plate.

Best wishes

Torsten.

[Walter Roberson]

1) suggest me to write the single function by combining both 'main' and 'D' functions.

Okay.

We suggest that you write a single function by combining both 'main' and 'D' functions, using parameters f, xa, xb, ya, and h

https://www.mathworks.com/help/matlab/ref/function.html

[PJS KUMAR]

I combined the two functions. But i am getting error. correct the code

function []=main4(fn,a,b,y0,h)

   x = a:h:b;
   y(1) = y0;

   for i = 1:numel(x)-1   
   dfstart= fn;
   dfold = dfstart;
   D(1) = double(subs(dfold,{x,y},{x(i),y(i)}));
   for i=2:n
     dfnew = subs(diff(dfold,x),diff(y(x),x),dfstart);
     D(i) =  double(subs(dfnew,{x,y},{x(i),y(i)}));
     dfold = dfnew;
   end
       y(i+1) = y(i)+h*D(1)+h^2/2*D(2)+h^3/6*D(3)+h^4/24*D(4)+h^5/120*D(5);
   end
  [x' y'] 
  end

I am getting the error as

>> main4(f,1,3,0.8,0.5)
Error using sym/subs>normalize (line 221)
Entries in second argument must be scalar.

Error in sym/subs>mupadsubs (line 147)
[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 135)
    G = mupadsubs(F,X,Y);

Error in main4 (line 10)
 D(1) = double(subs(dfold,{x,y},{x(i),y(i)}));

[Walter Roberson]

As I tried to get you to understand before:

If you are using a symbolic variable with a particular name, then you should strongly avoid using the same variable name for a different purpose, such as storing a list of particular locations to execute at, or such as storing the results of calculating the taylor series at particular locations.

[Torsten]

1. n is undefined.
2. You use i as a loop variable in both nested loops.
3. Walter's remark.

Answer 2

[Walter Roberson]

As I already explained to you, because you have

y = [ 3/2, x^2 + y(x)^2, 2*x + 2*y(x)*diff(y(x), x), 2*diff(y(x), x)^2 + 2*y(x)*diff(y(x), x, x) + 2]

then the y on the left side refers to the same thing as the y on the right side, and so y(x) on the right side signifies array indexing. Your y vector is 5 elements long, so the only valid values of x for this would be 1, 2, 3, 4, or 5. Each y(x) would resolve to a numeric scalar value, and diff() of a numeric scalar value is empty. Any operation involving an empty array returns an empty array, so all of the entries except the first two are going to disappear, leaving you only [3/2, x^2 + y(x)^2] to work with for preliminary consistency.

Now, with that subset, can x = 1 be made self-consistent? That would require that y(1) be 3/2, which seems plausible. With the subset, can x = 2 be made self-consistent? That would require that y(2) = x^2 + y(x)^2 = 2^2 + y(2)^2 . Rewriting as z = 4 + z^2 we can see that has only complex roots 1/2-1i*sqrt(15)*(1/2), 1/2+1i*sqrt(15)*(1/2) . Is that acceptable, to force y(2) to be complex valued?