Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
High Order Finite Difference

Subject: High Order Finite Difference

From: Lee

Date: 6 Mar, 2011 21:12:14

Message: 1 of 12

Hi,
Could someone please help.
I want to write a function m-file that when given an input argument for the Order of a Power Series Interpolating Polynomial, for example of the form:-

g(x) = a0*x^4 + a1*x^3 + a2*x^2 + a3*x + a4

this is for 5 data points, or, 4th Order Accurate.
Where the 5 data points are;
                     f0 f1 f2 f3 f4
                     X----------------X----------------X----------------X----------------X
                     x0 <---h---> x1 x2 x3 x4
                     0 h 2h 3h 4h
Will return an algebraic expression for the Finite Difference Approximations to the First and Second Derivatives.
For the above example will Return;
g'(x2) = g'(2h) = (f(i-2) - 8f(i-1) + 8f(i+1) - f(i+2))/(12h)
g''(x2) = g''(2h) = (- f(i-2) + 16f(i-1) - 30f(i) + 16f(i+1) - f(i+2))/(12h^2)
But these are the Central Difference Approximations, I need to be able to gain Expressions for all points, i.e. x1=h, x2=2h, x3=3h and x4=4h.
Also Not Just for a 4th Order Polynomial but Any Order Polynomial.

Every Program I have looked at so far only finds the coefficients of the polynomial, i.e. a0, a1, a2, a3, and a4. I know how to do this, but finding the Difference Approximations by hand is then very time consuming.
I need a program that will simply return the Approximations for any Order Polynomial for First and Second Derivatives and at Any Data Point.

If anyone has any ideas about how to do this could you please help...
also I am not too familiar with some of Matlabs more complex built in functions, so would prefer to use: for and while loops where possible.
I am also in abit of a hurry as I need the code for my final year Thesis, if someone could get back to me quickly it would be much appreciated.

Thank you
Lee

Subject: High Order Finite Difference

From: Nasser M. Abbasi

Date: 6 Mar, 2011 21:42:49

Message: 2 of 12

On 3/6/2011 1:12 PM, Lee wrote:
> Hi,

> I need a program that will simply return the Approximations for any Order
> Polynomial for First and Second Derivatives and at Any Data Point.
>

help taylor
help diff
help subs

I am sure with the above you can accomplish what you want.

--Nasser

Subject: High Order Finite Difference

From: Bruno Luong

Date: 6 Mar, 2011 22:20:19

Message: 3 of 12

"Lee " <369937@swansea.ac.uk> wrote in message <il0tbe$sru$1@fred.mathworks.com>...
> Hi,
> Could someone please help.
> I want to write a function m-file that when given an input argument for the Order of a Power Series Interpolating Polynomial, for example of the form:-
>
> g(x) = a0*x^4 + a1*x^3 + a2*x^2 + a3*x + a4
>
> this is for 5 data points, or, 4th Order Accurate.
> Where the 5 data points are;
> f0 f1 f2 f3 f4
> X----------------X----------------X----------------X----------------X
> x0 <---h---> x1 x2 x3 x4
> 0 h 2h 3h 4h
> Will return an algebraic expression for the Finite Difference Approximations to the First and Second Derivatives.
> For the above example will Return;
> g'(x2) = g'(2h) = (f(i-2) - 8f(i-1) + 8f(i+1) - f(i+2))/(12h)
> g''(x2) = g''(2h) = (- f(i-2) + 16f(i-1) - 30f(i) + 16f(i+1) - f(i+2))/(12h^2)
> But these are the Central Difference Approximations, I need to be able to gain Expressions for all points, i.e. x1=h, x2=2h, x3=3h and x4=4h.
> Also Not Just for a 4th Order Polynomial but Any Order Polynomial.
>
> Every Program I have looked at so far only finds the coefficients of the polynomial, i.e. a0, a1, a2, a3, and a4. I know how to do this, but finding the Difference Approximations by hand is then very time consuming.
> I need a program that will simply return the Approximations for any Order Polynomial for First and Second Derivatives and at Any Data Point.
>
> If anyone has any ideas about how to do this could you please help...
> also I am not too familiar with some of Matlabs more complex built in functions, so would prefer to use: for and while loops where possible.
> I am also in abit of a hurry as I need the code for my final year Thesis, if someone could get back to me quickly it would be much appreciated.
>

Try this:

p = 4; % polynomial order
h = 0.2; % step

x = (0:p)*h;
M = bsxfun(@power, x(:), (p:-1:0));

% point at which the derivative is to be calculated
x = x(3);
% first derivative coefficients
a = (p:-1:0) .* (x.^(p-1:-1:-1));
c1 = a/M;
% second derivative coefficients
a = (p:-1:0) .* (p-1:-1:-1) .* (x.^(p-2:-1:-2));
c2 = a/M;

% function values at x0-xp
f = rand(p+1,1);

fprintf('P''(x2) = %g\n', a1*f)
fprintf('P''''(x2) = %g\n', a2*f)

% Bruno

Subject: High Order Finite Difference

From: Lee

Date: 6 Mar, 2011 23:24:05

Message: 4 of 12

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <il11b3$ei5$1@fred.mathworks.com>...
> "Lee " <369937@swansea.ac.uk> wrote in message <il0tbe$sru$1@fred.mathworks.com>...
> > Hi,
> > Could someone please help.
> > I want to write a function m-file that when given an input argument for the Order of a Power Series Interpolating Polynomial, for example of the form:-
> >
> > g(x) = a0*x^4 + a1*x^3 + a2*x^2 + a3*x + a4
> >
> > this is for 5 data points, or, 4th Order Accurate.
> > Where the 5 data points are;
> > f0 f1 f2 f3 f4
> > X----------------X----------------X----------------X----------------X
> > x0 <---h---> x1 x2 x3 x4
> > 0 h 2h 3h 4h
> > Will return an algebraic expression for the Finite Difference Approximations to the First and Second Derivatives.
> > For the above example will Return;
> > g'(x2) = g'(2h) = (f(i-2) - 8f(i-1) + 8f(i+1) - f(i+2))/(12h)
> > g''(x2) = g''(2h) = (- f(i-2) + 16f(i-1) - 30f(i) + 16f(i+1) - f(i+2))/(12h^2)
> > But these are the Central Difference Approximations, I need to be able to gain Expressions for all points, i.e. x1=h, x2=2h, x3=3h and x4=4h.
> > Also Not Just for a 4th Order Polynomial but Any Order Polynomial.
> >
> > Every Program I have looked at so far only finds the coefficients of the polynomial, i.e. a0, a1, a2, a3, and a4. I know how to do this, but finding the Difference Approximations by hand is then very time consuming.
> > I need a program that will simply return the Approximations for any Order Polynomial for First and Second Derivatives and at Any Data Point.
> >
> > If anyone has any ideas about how to do this could you please help...
> > also I am not too familiar with some of Matlabs more complex built in functions, so would prefer to use: for and while loops where possible.
> > I am also in abit of a hurry as I need the code for my final year Thesis, if someone could get back to me quickly it would be much appreciated.
> >
>
> Try this:
>
> p = 4; % polynomial order
> h = 0.2; % step
>
> x = (0:p)*h;
> M = bsxfun(@power, x(:), (p:-1:0));
>
> % point at which the derivative is to be calculated
> x = x(3);
> % first derivative coefficients
> a = (p:-1:0) .* (x.^(p-1:-1:-1));
> c1 = a/M;
> % second derivative coefficients
> a = (p:-1:0) .* (p-1:-1:-1) .* (x.^(p-2:-1:-2));
> c2 = a/M;
>
> % function values at x0-xp
> f = rand(p+1,1);
>
> fprintf('P''(x2) = %g\n', a1*f)
> fprintf('P''''(x2) = %g\n', a2*f)
>
> % Bruno

Hi Bruno;
Cheers for the help..

I ran your code and got;

P'(x2)=1.90197
P''(x2)=4.80668

This is not what I wanted.
I need an algebraic expression.. for example. At Point: x2=2h, after running the code I want to see the following...

>> g'(x2) = g'(2h) = (f(i-2) - 8f(i-1) + 8f(i+1) - f(i+2))/(12h)
> > g''(x2) = g''(2h) = (- f(i-2) + 16f(i-1) - 30f(i) + 16f(i+1) - f(i+2))/(12h^2)

And similar expressions for any of the other points.
Thanks for trying though.

Lee

Subject: High Order Finite Difference

From: Lee

Date: 6 Mar, 2011 23:33:23

Message: 5 of 12

"Nasser M. Abbasi" <nma@12000.org> wrote in message <il0v4n$nor$1@speranza.aioe.org>...
> On 3/6/2011 1:12 PM, Lee wrote:
> > Hi,
>
> > I need a program that will simply return the Approximations for any Order
> > Polynomial for First and Second Derivatives and at Any Data Point.
> >
>
> help taylor
> help diff
> help subs
>
> I am sure with the above you can accomplish what you want.
>
> --Nasser

Hi Nasser;

I don't want to use a Taylor Series Expansion to derive the formulae, I need to use an Interpolating Polynomial, its kind of the whole idea of my Thesis.

Thanks anyway

Lee

Subject: High Order Finite Difference

From: Nasser M. Abbasi

Date: 6 Mar, 2011 23:49:43

Message: 6 of 12

On 3/6/2011 3:33 PM, Lee wrote:
> "Nasser M. Abbasi"<nma@12000.org> wrote in message<il0v4n$nor$1@speranza.aioe.org>...
>> On 3/6/2011 1:12 PM, Lee wrote:
>>> Hi,
>>
>>> I need a program that will simply return the Approximations for any Order
>>> Polynomial for First and Second Derivatives and at Any Data Point.
>>>
>>
>> help taylor
>> help diff
>> help subs
>>
>> I am sure with the above you can accomplish what you want.
>>
>> --Nasser
>

> Hi Nasser;
>
> I don't want to use a Taylor Series Expansion to derive the formulae, I need to use an
> Interpolating Polynomial, its kind of the whole idea of my Thesis.
>
> Thanks anyway
>
> Lee

Is the following demo something like what you want? it gives
finite difference scheme for any order for different number of points.

http://demonstrations.wolfram.com/FiniteDifferenceSchemesOfOneVariable/

But any way, if you can do it with paper and pencil, then I do not
see why you can not do it writing some symbolic manipulation code as well.

syms has all the commands to do what is needed.

good luck.

--Nasser

Subject: High Order Finite Difference

From: Bruno Luong

Date: 7 Mar, 2011 06:55:05

Message: 7 of 12

"Lee " <369937@swansea.ac.uk> wrote in message <il152l$ecv$1@fred.mathworks.com>...

>
> P'(x2)=1.90197
> P''(x2)=4.80668
>
> This is not what I wanted.
> I need an algebraic expression.. for example. At Point: x2=2h, after running the code I want to see the following...
>
> >> g'(x2) = g'(2h) = (f(i-2) - 8f(i-1) + 8f(i+1) - f(i+2))/(12h)
> > > g''(x2) = g''(2h) = (- f(i-2) + 16f(i-1) - 30f(i) + 16f(i+1) - f(i+2))/(12h^2)
>
> And similar expressions for any of the other points.
> Thanks for trying though.
>

Look closely how P'(x2) is computed: it's c1*f, meaning
c(1)*f(1) + c(2)*f(2) + .... c(p+1)*f(p+1)

This is the expression you want.

Bruno

Subject: High Order Finite Difference

From: Bruno Luong

Date: 7 Mar, 2011 07:02:06

Message: 8 of 12

"Nasser M. Abbasi" <nma@12000.org> wrote in message <il16io$9og$1@speranza.aioe.org>...

>
> But any way, if you can do it with paper and pencil, then I do not
> see why you can not do it writing some symbolic manipulation code as well.
>
> syms has all the commands to do what is needed.
>

Well very simple: that does not apply for people like me who do not have symbolic toolbox.

Bruno

Subject: High Order Finite Difference

From: Nasser M. Abbasi

Date: 7 Mar, 2011 07:07:28

Message: 9 of 12

On 3/6/2011 11:02 PM, Bruno Luong wrote:
> "Nasser M. Abbasi"<nma@12000.org> wrote in message<il16io$9og$1@speranza.aioe.org>...
>
>>
>> But any way, if you can do it with paper and pencil, then I do not
>> see why you can not do it writing some symbolic manipulation code as well.
>>
>> syms has all the commands to do what is needed.
>>
>

> Well very simple: that does not apply for people like me who do not have symbolic toolbox.
>
> Bruno

I see. I use student version, and Matlab student version always comes with
the basic symbolic toolbox build-in.

I guess commerical Matlab is different.

--Nasser

Subject: High Order Finite Difference

From: Bruno Luong

Date: 7 Mar, 2011 07:13:04

Message: 10 of 12

"Bruno Luong" <b.luong@fogale.findmycountry> wrote
>
> Look closely how P'(x2) is computed: it's c1*f, meaning
> c(1)*f(1) + c(2)*f(2) + .... c(p+1)*f(p+1)
>

To check with the expression you have provided here is the printout of the coefficients:

p = 4; % polynomial order
h = 0.2; % step

x = (0:p)*h;
M = bsxfun(@power, x(:), (p:-1:0));

% point at which the derivative is to be calculated
x = x(3);
% first derivative coefficients
a = (p:-1:0) .* (x.^(p-1:-1:-1));
c1 = a/M;
% second derivative coefficients
a = (p:-1:0) .* (p-1:-1:-1) .* (x.^(p-2:-1:-2));
c2 = a/M;

fprintf('\nc1 = %s\n', mat2str(c1,3));
yourcoefs1 = [1 -8 0 8 -1]/(12*h);
fprintf('yourcoefs1 = %s\n', mat2str(yourcoefs1,3));

fprintf('\nc2 = %s\n', mat2str(c2,3));
yourcoefs2 = [-1 16 -30 16 -1]/(12*h^2);
fprintf('yourcoefs2 = %s\n', mat2str(yourcoefs2,3));

% Bruno

Subject: High Order Finite Difference

From: Lee

Date: 7 Mar, 2011 12:22:04

Message: 11 of 12

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <il20i0$6ko$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote
> >
> > Look closely how P'(x2) is computed: it's c1*f, meaning
> > c(1)*f(1) + c(2)*f(2) + .... c(p+1)*f(p+1)
> >
>
> % Bruno

%=========================================
 
Cheers Bruno

It confused me because you had h=0.2; and I wanted h=1.
Essencially I want: h=h, i.e. stay symbolic.

Thanks alot, works great.. I still need to modify it for my purpose but is very useful.

Lee

Subject: High Order Finite Difference

From: Bruno Luong

Date: 7 Mar, 2011 12:46:21

Message: 12 of 12

>
> It confused me because you had h=0.2; and I wanted h=1.
> Essencially I want: h=h, i.e. stay symbolic.
>
> Thanks alot, works great.. I still need to modify it for my purpose but is very useful.

The first derivative (and coefficients) scaled with 1/h, the second derivative scaled with 1/h^2. So you can just leave h=1 in the calculation of c1 and c2, then scale later.

Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us