Replacing symbolic derivatives in MATLAB

Question

I have a symbolic function that looks like this

syms x y(x) h
fn(x) = y + (h^2*(diff(y(x), x) + 2))/2 + (h^5*diff(y(x), x, x, x, x))/120 + (h^3*diff(y(x), x, x))/6 + (h^4*diff(y(x), x, x, x))/24 + h*(2*x + y(x) - 1)

I'd like to replace all instances of derivatives of y with its first derivative, i.e.

subs(fn, sym('diff(y(x), x)'), dy)

where dy is already defined as

dy(x) = 2*x + y(x) - 1

The result is the following:

ans(x) =

y + (h^2*(2*x + y(x) + 1))/2 + (h^5*diff(y(x), x, x, x, x))/120 + (h^3*diff(y(x), x, x))/6 + (h^4*diff(y(x), x, x, x))/24 + h*(2*x + y(x) - 1)

It replaces the first derivative, but not the higher derivatives. What I want is for the h^5 term to have (h^5*diff(dy(x), x, x, x). Is there any way to do that?

My current method is pretty hackish and involves converting the sym to a string, replacing first derivatives with dy, then converting back to a sym and evaluating to reduce the order of each term of the series by one, but it has to be recursive because at each stage derivatives of dy are then replaced by something containing diff(y, ...). I was hoping there would be a cleaner way to deal with this.

Answer 1

You need to keep in mind that Matlab treats things like diff(y,x) and diff(y,x,2) as distinct variables. It doesn't know how to substitute diff(y,x) into diff(y,x,2) because such a general operation for an abstract function (one with out an explicit definition, e.g., y(x)) is ill-defined.

How about something like this that performs substitution from the opposite end, starting with the highest order derivatives:

syms y(x) h
dy(x) = 2*x + y - 1
fn(x) = y + (h^2*(diff(y, x) + 2))/2 + (h^5*diff(y, x, 4))/120 + (h^3*diff(y, x, 2))/6 + (h^4*diff(y, x, 3))/24 + h*(2*x + y - 1)
fn2 = subs(fn, diff(y, x, 4), diff(dy, x, 3));
fn2 = subs(fn2, diff(y, x, 3), diff(dy, x, 2));
fn2 = subs(fn2, diff(y, x, 2), diff(dy, x));
fn2 = subs(fn2, diff(y, x), dy);

This returns

fn2(x) =

y(x) + (h^2*(2*x + y(x) + 1))/2 + (h^3*(2*x + y(x) + 1))/6 + (h^4*(2*x + y(x) + 1))/24 + (h^5*(2*x + y(x) + 1))/120 + h*(2*x + y(x) - 1)

Or you can leave dy(x) as an abstract symbolic expression initially:

syms y(x) dy(x) h
fn(x) = y + (h^2*(diff(y, x) + 2))/2 + (h^5*diff(y, x, 4))/120 + (h^3*diff(y, x, 2))/6 + (h^4*diff(y, x, 3))/24 + h*(2*x + y - 1)
fn2 = subs(fn, diff(y, x, 4), diff(dy, x, 3));
fn2 = subs(fn2, diff(y, x, 3), diff(dy, x, 2));
fn2 = subs(fn2, diff(y, x, 2), diff(dy, x));
fn2 = subs(fn2, diff(y, x), dy)

which returns

fn2(x) =

y(x) + (h^4*diff(dy(x), x, x))/24 + (h^2*(dy(x) + 2))/2 + (h^5*diff(dy(x), x, x, x))/120 + (h^3*diff(dy(x), x))/6 + h*(2*x + y(x) - 1)