Calculate derivatives on the interval 0<x<2pi using forward, backward, centered, and fourth order accurate finite differences

f(x) = ax+bx^2+cx^3+dsin(mx)+gexp(hx)

a=-5
b=.5
c=-.15
d=.2
m=6
g=-2
h=-3

Calculate first derivatives using forward, backward, centered, and fourth order accurate finite differences.

calculate second derivative using 2nd order accurate central differences.

plot the finite difference solutions compared to the analytic solutions. 6 plots total (3 for the first derivatives and 3 for the second derivatives) with 5 lines per graph for the first derivatives and 2 lines per graph for the second derivatives.

for simplicity it is ok to plot the solutions starting at grid point 3, and ending at grid point N-2, or any alternate methods for calculating the derivatives near the boundaries may be used that is consistent with the interior scheme.

also calculate the RMS error for each finite difference derivative at each grid resolution and present the results in a table or graph showing N vs error for each scheme. Note that these should form nearly a straight line on a log-log scale with a slope of -1 for the first order accurate schemes and a slope of -2 for the second order accurate schemes and a slope of -4 for the fourth order accurate schemes.

This is what I've done so far...but I don't know where to go from here or what I'm doing wrong.

    clear
  hw4_input
  L = 2*pi;
  N=input('How many grid points would you like to use on the interval. Enter 20,200, or 2000');
  if (N==20);
      dx = L/(N-1);
      x=0:dx:2*pi;
      analytical
          for i = 1:N
              f(i) = a*x(i)+b*x(i)^2+c*x(i)^3+d*sin(m*x(i))+g*exp(h*x(i));
              dfdx=zeros(1,N);
              dfdx(i) = a+2*b*x(i)+3*c*x(i)^2+m*d*cos(m*x(i))+h*g*exp(h*x(i));
              x(i)=x(i);
          end
          fw
          for i = 1:N-1
              f(i) = a*x(i)+b*x(i)^2+c*x(i)^3+d*sin(m*x(i))+g*exp(h*x(i));
              dfdx1=zeros(1,N);
              dfdx1(i) = (f(i+1)-f(i))/dx;
              x1(i) = x(i);
          end
          bw
          for i= 2:N
              f(i) = a*x(i)+b*x(i)^2+c*x(i)^3+d*sin(m*x(i))+g*exp(h*x(i));
              dfdx2=zeros(1,N);
              dfdx2(i) = (f(i)-f(i-1))/dx;
              x2(i) = x(i);
          end
          cd
          for i = 2:N-1
              f(i) = a*x(i)+b*x(i)^2+c*x(i)^3+d*sin(m*x(i))+g*exp(h*x(i));
              dfdx3=zeros(1,N);
              dfdx3(i)=(f(i+1)-f(i-1))/(2*dx);
              x3(i)= x(i);
          end
          forth
          for i = 1:N-3
              f(i) = a*x(i)+b*x(i)^2+c*x(i)^3+d*sin(m*x(i))+g*exp(h*x(i));
              dfdx4=zeros(N,1);
              dfdx4(i)= (-f(i+2)+8*f(i+1)-8*f(i-1)+f(i-2))/(12*dx);
          end
  end

Answer 1

[Matt Tearle]

First, those for loops are making Baby Cleve cry. They can all be replaced with vectorized expressions ( f = a*x + b*x.^2 + etc).

Second, I don't see where you've defined the constants a, b, c, etc EDIT: is that what hw4_input does? In which case, never mind.

Third, your comments (analytical, fw, bw, etc) don't seem to commented, but maybe that's just a cut-paste issue.

Fourth, I don't see why you have the if statement - it should work for any N

Fifth, differences can be calculated in a vectorized manner using diff or indexing (and vector subtraction)

Sixth, plot (and related log scale plot functions) allow multiple x-y pairs. You'll need to use indexing to extract the appropriate elements of x. The colon operator, to define index ranges, is your friend - eg x(4:end) extracts the fourth through the last elements of x.

EDIT: Seventh, your commands dfdx=zeros(1,N); inside the loops overwrite everything on each iteration! As a result, you'll be left with a vector of zeros except for one value (the last run through the loop). If you get rid of the loops, it won't matter, but if you feel you must keep the loops, at least move these preallocation commands before the for statements.

Let's see if that gets you any further forward...

[Matt Tearle]

I have more, but let's get this much working first...

When is this due...?

Answer 2

[Matt Fig]

It looks like you just pasted the HW assignment. But this is MATLAB Answers, implying that there are questions being answered. I don't see any questions in your paste.

Soooo, do you have a question?

[Austin]

haha yeah I put up what I've done so far but don't know where to go from there or how to plot it all