Taylor series method expansion

function sol=taylor1(fn,dfn,a,b,y0,n)
h=(b-a)/n;
x=a+(0:n)*h;
y(1)=y0;
for k=1:n
    y(k+1)=y(k)+h*feval(fn,x(k),y(k)) + (h^2*feval(dfn,x(k),y(k)))/2;
end
sol=[x',y'];

the above function working for the functions with single variable 'x', how can we evaluate for two variables x and y. Giving error as

>> taylor1(f,df,0,0.8,1.5,5)
Error using feval
Argument must contain a character vector or function handle.
Error in taylor1 (line 8)
    y(k+1)=y(k)+h*feval(fn,x(k),y(k)) + (h^2*feval(dfn,x(k),y(k)))/2;

2) how can we expand the series for more number of terms 3) how can we use vectors instead of for loop

[Walter Roberson]

What is the f that you are passing as the first argument to taylor1 ?

[PJS KUMAR]

derivative of y in terms of x and y

[Walter Roberson]

What form is the f that you are passing in? Is it a function handle? Is a symbolic expression? Is it a double?

I suspect that you are passing in a symbolic expression. You cannot feval() a symbolic expression: you need to subs() into it. Or you need to use matlabFunction() to convert the symbolic expression into an anonymous function.

[PJS KUMAR]

I am passing a function handle

[Walter Roberson]

Please show your code. Because by the time of the feval, what has reached fn is not a function handle.

[PJS KUMAR]

see the code. error in line 8

[Walter Roberson]

That is not the same error message as before.

The problem now is that your invocation of feval() passes x(k) and y(k) to the provided function, but your provided function expects only a single parameter. If you wish to ignore the second parameter you can do that but you still have to code it, such as

fn = @(x,y) x^3+1

[Titus Edelhofer]

The problem is not line 8 but how you call your taylor1. You have something like

f = @(x,y) something;
df = @(x,y) somethingelse;

and then you call taylor1. The error is in the df=... line. But we can't help since we haven't seen it ...

Titus

[PJS KUMAR]

suggest the code for the following taylor series expansion, given initial values x0, f(x0)

1. one vector should contain the values of f(x0), f'(x0), f''(x0), f'''(x0),....
2. second vector should contain the values of 1/0!, (x-x0)/1!, (x-x0)^2/2!, (x-x0) ^3/3!,....

[PJS KUMAR]

Suggest me how to write the following code using vectors without using for loop

    y(1)=y0;
    for k=1:n
        y(k+1)=y(k)+h*feval(fn,x(k),y(k));
    end

[PJS KUMAR]

given x0,y0,h,x values

can we write the following code using x and y vectors

    x1=x0;
    y1=y0;
    while(x>x1)
      y1=y1+h.*df(x1,y1);
      x1=x1+h;
      fprintf('When x=%2.2f y=%2.2f\n',x1,y1)
    end

[PJS KUMAR]

I have the following terms

y'=x-y^2
y''= diff(y',x,1)
y'''= diff(y',x,2)
y''''=diff(y',x,3)

and so on

can we represent these values in the form of a vector

[Walter Roberson]

"one vector should contain the values of f(x0), f'(x0), f''(x0), f'''(x0),...."

"second vector should contain the values of 1/0!, (x-x0)/1!, (x-x0)^2/2!, (x-x0) ^3/3!,"

There are some non-loop possibilities from there. Assuming they are row vectors and not column vectors:

• sum(first_vector .* second_vector)
• first_vector * second_vector.'

In the case of real valued functions, you can also do

• dot(first_vector, second_vector)

Answer 1

[PJS KUMAR]

I wrote the following code for tayler 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:5)./factorial(0:5)
  for i=2:5
        y(i+1)=diff(y(i),x)
  end
  for i=1:5
  x=a+i*h;
  c=eval(y);
  yb=sum(c.*ht);
  fprintf('The value of %f is %f\n',x,yb);
  end
  end

1) can we change the code for n order (n derivatives) as i am using upto 5

2) suggest me to change the code for the functions with two variables

3) can we use vectorisation to avoid 'for' loops

[Walter Roberson]

Order n:

function yb=taylor(f,a,b,ya,n)
  syms x y;
  h=(b-a)/n;
  y(1)=ya;
  y(2)=f;
  ht=h.^(0:n)./factorial(0:n)
  for i=2:n
        y(i+1)=diff(y(i),x)
  end
  for i=1:n
  x=a+i*h;
  c=eval(y);
  yb=sum(c.*ht);
  fprintf('The value of %f is %f\n',x,yb);
  end
  end

2) The code already handles functions in multiple variables, as long as the variable of differentiation is x. To change that you would have to pass in the variable of differentiation and use that instead of x in your diff() call. Notice, though, that whatever your a, b, and ya parameters are supposed to mean, they have been constructed to only deal with a single dimension, not to deal with the possibility of expansion around a 2D point or in a circle or square.

Reminder: taylor series are always only with respect to one variable at a time. In some contexts it can make sense to extend taylor series to multiple dimensions around a point. In such case the way to proceed is to take the taylor series with respect to a single variable at a time.

3) There is no way to ask for differentiation to a number of different degrees simultaneously, only to ask for differentiation to one particular degree. It is therefor not possible to vectorize the differentiation process. The closest you could get would be to use arrayfun() to hide the differentiation loop, which would probably involve recalculating the same derivative many times unless you took care to "memoize" the calculation.

4) Your code will fail no matter what the data type is that you give for ya and f. The closest you could get would be if you specified both ya and f as symbolic, then you could get as far as the eval(). However, eval of a symbolic array can fail, because it means the same thing as eval() of char() of the symbolic array except replacing matrix() and array() calls with better syntax, and char() of symbolic entries can involve syntax that is not MATLAB and is not MuPAD either. You should never eval() a symbolic expression.

Answer 2

[PJS KUMAR]

The above program is not working for the function f=x^2 + y^2 and giving error as

>> f=@(x)x^2+y^2
f =
  function_handle with value:
    @(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)]
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]
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, 2*y(x)*diff(y(x), x, x, x) + 6*diff(y(x), x)*diff(y(x), x, x)]
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, 2*y(x)*diff(y(x), x, x, x) + 6*diff(y(x), x)*diff(y(x), x, x), 2*y(x)*diff(y(x), x, x, x, x) + 6*diff(y(x), x, x)^2 + 8*diff(y(x), x)*diff(y(x), x, x, x)]
Subscript indices must either be real positive integers or logicals.
Error in sym/eval (line 11)
s = evalin('caller',vectorize(map2mat(char(x))));
Error in taylor (line 12)
c=eval(y)

Suggest necessary corrections in the code. Also tell me how to evaluate y for all values of x in between a and b

[PJS KUMAR]

I am calling the function taylor(f,a,b,ya,n) with f=x^2+y^2, a=0, b=0.8, ya=1.5, n=4

>> f=@(x)x^2+y^2
f =
  function_handle with value:
    @(x)x^2+y^2
>> taylor(f,0,0.8,1.5,4)

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:5
      y(i+1)=diff(y(i),x)
end

and created two vectors 'ht' and 'y' where

ht=(1,h,h^2/2!,h^3/3!,....)

y=(ya,y',y'',y''',........)

and obtained the result as

• 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, 2*y(x)*diff(y(x), x, x, x) + 6*diff(y(x), x)*diff(y(x), x, x), 2*y(x)*diff(y(x), x, x, x, x) + 6*diff(y(x), x, x)^2 + 8*diff(y(x), x)*diff(y(x), x, x, x)]

give me the code to complete the function with the following actions

• I want to evaluate y at different values of x ranging from a to b, i.e., to evaluate y',y'',y'''.... at values of x and to store the values in vector c.
• To multiply the elements of 'ht' and 'c' vectors and to display the sum of the products

[Walter Roberson]

You have a conflict between the variable named y in your code in the form

 y(1)=ya;
 y(2)=f;

and the y in your f=@(x)x^2+y^2, which refers to some variable that you must have defined before you called the function. The conflict exists even if you had used

syms y
f = @(x)x^2+y^2

After all, look at what you are saying: 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, 2*y(x)*diff(y(x), x, x, x) + 6*diff(y(x), x)*diff(y(x), x, x), 2*y(x)*diff(y(x), x, x, x, x) + 6*diff(y(x), x, x)^2 + 8*diff(y(x), x)*diff(y(x), x, x, x)]

When you eval() this, y is the vector shown, so y(x) is asking to index the vector at the location given by x, which you have defined as x=a+i*h which is unlikely to work out as an exact integer.

If you had used some other variable name instead of y(1) and so on, then y in the eval() would instead refer to

syms x y;

which would be a scalar, and you would have the related problem of trying to index a scalar symbol at a non-integer location (in particular at a location that did not happen to exactly equal 1)

We can tell from the output you show that you must have done

syms y(x)

before you defined

f = @(x)x^2+y^2

But is it really true that y is a function of x, that the expression to be taylor'd is x^2 + y(x)^2 ?

Anyhow, diff(y(x), x, x) when x is numeric and y(x) is symbolic is not defined. In such a case, diff(y(x), x) when x is numeric and y(x) is symbolic is defined, but only for the case where x happens to be a non-negative integer, in which case it is the derivative number. And if y(x) is numeric instead of symbolic and x is numeric, then diff(y(x), x, x) is only valid if x is a nonnegative integer dimension number (third parameter) and is simultaneously a nonnegative integer designating the order difference (second parameter).

Remember I posted above that you should never eval() a symbolic expression? This is an example of what happens: you end up mixing symbolic expression syntax with MATLAB syntax, and you get Fail City.

[Elavarasi E]

am too working on a similar kind of problem. can u guide me if ur code has worked