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:
curvature radius

Subject: curvature radius

From: Junghun

Date: 20 Oct, 2010 02:21:03

Message: 1 of 9

Hi,
I've been trying to find the curvature radius of my data set.
x1=[x1(1) x1(2)....]
y1=[y1(1) y1(2)....]

x1 is just numbers, and the interval between x1(i) and x1(i+1) is constant
like x1= [1 3 5 7...]. dx1=2

y1 is a vector (y coordinate of x1) which are variables.

what I did (for getting curvature radius for each point) is as follows.
suppose we wanna get curvature at x1(2)

k1=(y1(3)-y1(1))/dx1 -> 1st derivative
k2=((y1(3)-y1(2))/dx1-(y1(2)-y1(1))/dx1)/dx1 -> 2nd derivative

r= (1+k2^2)^(3/2)/abs(k2)

In this way, I could get r in term of y1(1), y1(2), and y1(3).

My problem is as follows.
If I change dx1 value, then I will get totally different function for the curvature radius r.
I don't know how to resolve this scaling problem.

Is there any way to resolve this problem?

p.s.
When I used one of the other methods for getting 'r' with using solving the least square problem, it didn't calculate r correctly from the first place. I don't know why...


  

Subject: curvature radius

From: Roger Stafford

Date: 20 Oct, 2010 03:18:04

Message: 2 of 9

"Junghun " <jcho7@wisc.edu> wrote in message <i9ljmf$1gm$1@fred.mathworks.com>...
> Hi,
> I've been trying to find the curvature radius of my data set.
> x1=[x1(1) x1(2)....]
> y1=[y1(1) y1(2)....]
>
> x1 is just numbers, and the interval between x1(i) and x1(i+1) is constant
> like x1= [1 3 5 7...]. dx1=2
>
> y1 is a vector (y coordinate of x1) which are variables.
>
> what I did (for getting curvature radius for each point) is as follows.
> suppose we wanna get curvature at x1(2)
>
> k1=(y1(3)-y1(1))/dx1 -> 1st derivative
> k2=((y1(3)-y1(2))/dx1-(y1(2)-y1(1))/dx1)/dx1 -> 2nd derivative
>
> r= (1+k2^2)^(3/2)/abs(k2)
>
> In this way, I could get r in term of y1(1), y1(2), and y1(3).
>
> My problem is as follows.
> If I change dx1 value, then I will get totally different function for the curvature radius r.
> I don't know how to resolve this scaling problem.
>
> Is there any way to resolve this problem?
>
> p.s.
> When I used one of the other methods for getting 'r' with using solving the least square problem, it didn't calculate r correctly from the first place. I don't know why...
- - - - - - - - - -
  You should not be surprised that rescaling x without rescaling y by the same factor will make changes on the radius which are disproportionate. If you perform the same rescaling on both x and y values, the radius would also be rescaled in the same way, but if they are rescaled differently the center of curvature will be displaced and the radius altered in a more complicated way than could be accounted for by any kind of rescaling.

  I would think you would be happier drawing the unique circle through the three points to determine the radius of curvature. You can find the center of curvature at the intersection of the two perpendicular bisectors of segments P1P2 and P2P3 where P1 = (x1(1),y1(1)), P2 = (x1(2),y1(2)), and P3 = (x1(3),y1(3)). From there you can find the radius. Matlab is good at solving such a problem.

  I see a couple of errors in your formulas. For k1 you should have

k1=(y1(3)-y1(1))/(2*dx1)

For r you should have

r= (1+k1^2)^(3/2)/abs(k2)

Roger Stafford

Subject: curvature radius

From: Junghun

Date: 20 Oct, 2010 03:41:03

Message: 3 of 9

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9ln1c$s6h$1@fred.mathworks.com>...
> "Junghun " <jcho7@wisc.edu> wrote in message <i9ljmf$1gm$1@fred.mathworks.com>...
> > Hi,
> > I've been trying to find the curvature radius of my data set.
> > x1=[x1(1) x1(2)....]
> > y1=[y1(1) y1(2)....]
> >
> > x1 is just numbers, and the interval between x1(i) and x1(i+1) is constant
> > like x1= [1 3 5 7...]. dx1=2
> >
> > y1 is a vector (y coordinate of x1) which are variables.
> >
> > what I did (for getting curvature radius for each point) is as follows.
> > suppose we wanna get curvature at x1(2)
> >
> > k1=(y1(3)-y1(1))/dx1 -> 1st derivative
> > k2=((y1(3)-y1(2))/dx1-(y1(2)-y1(1))/dx1)/dx1 -> 2nd derivative
> >
> > r= (1+k2^2)^(3/2)/abs(k2)
> >
> > In this way, I could get r in term of y1(1), y1(2), and y1(3).
> >
> > My problem is as follows.
> > If I change dx1 value, then I will get totally different function for the curvature radius r.
> > I don't know how to resolve this scaling problem.
> >
> > Is there any way to resolve this problem?
> >
> > p.s.
> > When I used one of the other methods for getting 'r' with using solving the least square problem, it didn't calculate r correctly from the first place. I don't know why...
> - - - - - - - - - -
> You should not be surprised that rescaling x without rescaling y by the same factor will make changes on the radius which are disproportionate. If you perform the same rescaling on both x and y values, the radius would also be rescaled in the same way, but if they are rescaled differently the center of curvature will be displaced and the radius altered in a more complicated way than could be accounted for by any kind of rescaling.
>
> I would think you would be happier drawing the unique circle through the three points to determine the radius of curvature. You can find the center of curvature at the intersection of the two perpendicular bisectors of segments P1P2 and P2P3 where P1 = (x1(1),y1(1)), P2 = (x1(2),y1(2)), and P3 = (x1(3),y1(3)). From there you can find the radius. Matlab is good at solving such a problem.
>
> I see a couple of errors in your formulas. For k1 you should have
>
> k1=(y1(3)-y1(1))/(2*dx1)
>
> For r you should have
>
> r= (1+k1^2)^(3/2)/abs(k2)
>
> Roger Stafford


Thanks, Roger.

Could you be more specific about the method you mentioned?
Is the method same as what you mentioned in "curvature and radius of curvature of a plane" by using the code below?

mx = mean(x); my = mean(y)
 X = x - mx; Y = y - my; % Get differences from means
 dx2 = mean(X.^2); dy2 = mean(Y.^2); % Get variances
 t = [X,Y]\(X.^2-dx2+Y.^2-dy2)/2; % Solve least mean squares problem
 a0 = t(1); b0 = t(2); % t is the 2 x 1 solution array [a0;b0]
 r = sqrt(dx2+dy2+a0^2+b0^2); % Calculate the radius
 a = a0 + mx; b = b0 + my; % Locate the circle's center
 curv = 1/r; % Get the curvature

I used this code, but it didn't work. I don't the exact reason. I just guess it's because the r is including variables. In my actual code, the follows is my "flow".

1) Get r in term of y variables
2) Get f(r)
3) find the set of y, which makes f(r) minimize.

Or, the method that you mentioned is different from what we mentioned in "curvature and radius of curvature of a plane"?

Subject: curvature radius

From: Roger Stafford

Date: 20 Oct, 2010 07:54:03

Message: 4 of 9

"Junghun " <jcho7@wisc.edu> wrote in message <i9locf$mlf$1@fred.mathworks.com>...
> Thanks, Roger.
>
> Could you be more specific about the method you mentioned?
> Is the method same as what you mentioned in "curvature and radius of curvature of a plane" by using the code below?
>
> mx = mean(x); my = mean(y)
> X = x - mx; Y = y - my; % Get differences from means
> dx2 = mean(X.^2); dy2 = mean(Y.^2); % Get variances
> t = [X,Y]\(X.^2-dx2+Y.^2-dy2)/2; % Solve least mean squares problem
> a0 = t(1); b0 = t(2); % t is the 2 x 1 solution array [a0;b0]
> r = sqrt(dx2+dy2+a0^2+b0^2); % Calculate the radius
> a = a0 + mx; b = b0 + my; % Locate the circle's center
> curv = 1/r; % Get the curvature
>
> I used this code, but it didn't work. I don't the exact reason. I just guess it's because the r is including variables. In my actual code, the follows is my "flow".
>
> 1) Get r in term of y variables
> 2) Get f(r)
> 3) find the set of y, which makes f(r) minimize.
>
> Or, the method that you mentioned is different from what we mentioned in "curvature and radius of curvature of a plane"?
- - - - - - - - - -
  The method you referred to that I wrote in July of 2007 in the thread "curvature and radius of curvature of a plane" was actually intended for finding a least squares fit to a circle for a large number of points. However when applied to just three points it does give the center and radius of the precise circle that runs through those three points. I have just successfully tested it this evening for a number of randomly placed triplets of points.

  However it is a rather cumbersome method of finding the circle for only three points. Below is a somewhat more direct way of solving the three-point problem. Let the three points be (x1,y1), (x2,y2), and (x3,y3).

t1 = x3^2-x2^2+y3^2-y2^2;
t2 = x1^2-x3^2+y1^2-y3^2;
t3 = x2^2-x1^2+y2^2-y1^2;
d = x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3;
a = 1/2*(t1*y1+t2*y2+t3*y3)/d;
b = -1/2*(t1*x1+t2*x2+t3*x3)/d;

Then (a,b) will be the center of a circle running through the three points. I have tested this one too and it gives the same answer as the first method. In both methods the distances from (a,b) to the three points are all equal (allowing for round off error.)

  You say about the first method that "it didn't work". I challenge you to try both of these methods on some point triplets and see if they don't actually find points which are equidistant from the three points in every case. If you find an example that fails, please show me the three points. Of course the three points must not be colinear, since their curvature would then be zero and no solution would be possible.

Roger Stafford

Subject: curvature radius

From: Roger Stafford

Date: 20 Oct, 2010 14:28:05

Message: 5 of 9

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9m76r$nr7$1@fred.mathworks.com>...
> .......
> However it is a rather cumbersome method of finding the circle for only three points. Below is a somewhat more direct way of solving the three-point problem. Let the three points be (x1,y1), (x2,y2), and (x3,y3).
>
> t1 = x3^2-x2^2+y3^2-y2^2;
> t2 = x1^2-x3^2+y1^2-y3^2;
> t3 = x2^2-x1^2+y2^2-y1^2;
> d = x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3;
> a = 1/2*(t1*y1+t2*y2+t3*y3)/d;
> b = -1/2*(t1*x1+t2*x2+t3*x3)/d;
>
> Then (a,b) will be the center of a circle running through the three points. I have tested this one too and it gives the same answer as the first method. In both methods the distances from (a,b) to the three points are all equal (allowing for round off error.)
> ........
- - - - - - - - - - -
  I should have mentioned that finding the circle through the three points is simply the problem of finding the circumscribed circle about a triangle with the three points as vertices. The radius of such a circle by a well-known formula is the product of the three sides divided by four times the area of the triangle.

Roger Stafford

Subject: curvature radius

From: Junghun

Date: 21 Oct, 2010 04:17:11

Message: 6 of 9

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9mu9l$alk$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9m76r$nr7$1@fred.mathworks.com>...
> > .......
> > However it is a rather cumbersome method of finding the circle for only three points. Below is a somewhat more direct way of solving the three-point problem. Let the three points be (x1,y1), (x2,y2), and (x3,y3).
> >
> > t1 = x3^2-x2^2+y3^2-y2^2;
> > t2 = x1^2-x3^2+y1^2-y3^2;
> > t3 = x2^2-x1^2+y2^2-y1^2;
> > d = x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3;
> > a = 1/2*(t1*y1+t2*y2+t3*y3)/d;
> > b = -1/2*(t1*x1+t2*x2+t3*x3)/d;
> >
> > Then (a,b) will be the center of a circle running through the three points. I have tested this one too and it gives the same answer as the first method. In both methods the distances from (a,b) to the three points are all equal (allowing for round off error.)
> > ........
> - - - - - - - - - - -
> I should have mentioned that finding the circle through the three points is simply the problem of finding the circumscribed circle about a triangle with the three points as vertices. The radius of such a circle by a well-known formula is the product of the three sides divided by four times the area of the triangle.
>
> Roger Stafford

Thanks, Roger
Here I attached the code below.

In this case, r is not what we expected in the sense that:
1) we cannot see the first peak at around 0.2*10^-3
2) the height of peaks is different (getting bigger) which should be equal in principle.
3) If you change dx, dx1 you will see different pattern, which means r is sensitive to dx and dx1.

I know this is sine curve, which has zero slope, which makes troubles.
If so, Is there any way to resolve this kind of problem?

The code===================================

clear all;

x=(4*pi)/(20*10.^4);
dx=(4*pi)/(20*10.^4);
dx1=(4*pi)/(20*10.^4);

for i=1:21
    y0=5*10.^-8*sin(10.^4*x);
    y1(i)=y0;
    x1(i)=x;
    x=x+dx;
    
end


for j=1:21
   if (j==1)
        k1(2)=y1(20);
        k2(2)=y1(j);
        k3(2)=y1(j+1);
        
        k1(1)=x1(j)-dx1;
        k2(1)=x1(j);
        k3(1)=x1(j+1);
        
    else if (j==21)
            k1(2)=y1(j-1);
            k2(2)=y1(j);
            k3(2)=y1(2);
            
            k1(1)=x1(j-1);
            k2(1)=x1(j);
            k3(1)=x1(j)+dx1;
            
        else
           
            k1(2)=y1(j-1);
            k2(2)=y1(j);
            k3(2)=y1(j+1);
            
            k1(1)=x1(j-1);
            k1(1)=x1(j);
            k1(1)=x1(j+1);
            
        end
   end
        k1(3)=0;
        k2(3)=0;
        k3(3)=0;
        
        % using "Fit circle to 3 points" method
        tt=k2-k1;
        uu=k3-k1;
        vv=k3-k2;
        ww=cross(tt,uu);
        tt2=sum(tt.^2);
        uu2=sum(uu.^2);
        vv2=sum(vv.^2);
        ww2=sum(ww.^2);
        
        r=1/2*sqrt(tt2*uu2*vv2/ww2);
               
        r1(j)=r;
   
end

r1
figure(1); clf;
plot(x1,r1,'o-')
figure(2); clf;
plot(x1,y1,'o-')

Subject: curvature radius

From: Roger Stafford

Date: 21 Oct, 2010 21:56:04

Message: 7 of 9

"Junghun " <jcho7@wisc.edu> wrote in message <i9oes7$na4$1@fred.mathworks.com>...
> .........
> k1(1)=x1(j-1);
> k1(1)=x1(j);
> k1(1)=x1(j+1);
> .........
- - - - - - - - - - -
  If the above is the code you actually used, you might have seen some strange results. Notice that for 1 < j < 21, k2(1) and k3(1) remain unchanged from the values they had at j = 1. I am sure this is not what you intended. Didn't you mean this:?

            k1(1)=x1(j-1);
            k2(1)=x1(j);
            k3(1)=x1(j+1);

  Just for the "record", for finding both a circle's center and radius the following requires fewer operations than the code I gave earlier. Again the three given points on the circle are (x1,y1), (x2,y2), and (x3,y3).

x21 = x2-x1; y21 = y2-y1;
x31 = x3-x1; y31 = y3-y1;
h21 = x21^2+y21^2; h31 = x31^2+y31^2;
d = 2*(x21*y31-x31*y21);
a = x1+(h21*y31-h31*y21)/d;
b = y1-(h21*x31-h31*x21)/d;
r = sqrt(h21*h31*((x3-x2)^2+(y3-y2)^2))/abs(d);

The circle's center is at (a,b) and its radius is r.

Roger Stafford

Subject: curvature radius

From: Junghun

Date: 15 Nov, 2010 04:21:03

Message: 8 of 9

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9qctk$cat$1@fred.mathworks.com>...
> "Junghun " <jcho7@wisc.edu> wrote in message <i9oes7$na4$1@fred.mathworks.com>...
> > .........
> > k1(1)=x1(j-1);
> > k1(1)=x1(j);
> > k1(1)=x1(j+1);
> > .........
> - - - - - - - - - - -
> If the above is the code you actually used, you might have seen some strange results. Notice that for 1 < j < 21, k2(1) and k3(1) remain unchanged from the values they had at j = 1. I am sure this is not what you intended. Didn't you mean this:?
>
> k1(1)=x1(j-1);
> k2(1)=x1(j);
> k3(1)=x1(j+1);
>
> Just for the "record", for finding both a circle's center and radius the following requires fewer operations than the code I gave earlier. Again the three given points on the circle are (x1,y1), (x2,y2), and (x3,y3).
>
> x21 = x2-x1; y21 = y2-y1;
> x31 = x3-x1; y31 = y3-y1;
> h21 = x21^2+y21^2; h31 = x31^2+y31^2;
> d = 2*(x21*y31-x31*y21);
> a = x1+(h21*y31-h31*y21)/d;
> b = y1-(h21*x31-h31*x21)/d;
> r = sqrt(h21*h31*((x3-x2)^2+(y3-y2)^2))/abs(d);
>
> The circle's center is at (a,b) and its radius is r.
>
> Roger Stafford

Thanks Roger,
I am having same kind of problems that I asked before.

Now, a curve is given. say, y=5.05*10.^-8*sin(10.^4*x)+0.4302*10.^-6
I hope to get same curvature radius, r graph even if I change dx scale.

I attach my code for that. (which is rough.)
In the code below, I set 4 different x intervals, and got the r graph.
The curvature radius is denpending on dx. (which I don't want to get)

Even if I tried the vector calculus that you taught me,
(by finding the radius of a circle which includes 3 given points, with vector calculus)
 r is still depending on dx.

In the fixed region, I hope to get the same r graph even when dx is different.
Could you help me out with this?



=======================================
clear all;

x1 =(6*pi)/(40*10.^4):(6*pi)/(40*10.^4):(6*pi)/(40*10.^4)*41;
x2 =(6*pi)/(30*10.^4):(6*pi)/(30*10.^4):(6*pi)/(30*10.^4)*31;
x3 =(6*pi)/(20*10.^4):(6*pi)/(20*10.^4):(6*pi)/(20*10.^4)*21;
x4 =(6*pi)/(10000*10.^4):(6*pi)/(10000*10.^4):(6*pi)/(10000*10.^4)*10001;

y1 = 5.05*10.^-8*sin(10.^4*x1)+0.4302*10.^-6;
y2 = 5.05*10.^-8*sin(10.^4*x2)+0.4302*10.^-6;
y3 = 5.05*10.^-8*sin(10.^4*x3)+0.4302*10.^-6;
y4 = 5.05*10.^-8*sin(10.^4*x4)+0.4302*10.^-6;

dx1 =(6*pi)/(40*10.^4)
dx2 =(6*pi)/(30*10.^4)
dx3 =(6*pi)/(20*10.^4)
dx4 =(6*pi)/(10000*10.^4)

for j=1:41 % both end points ignored
   if (j==1)
        k1=y1(40);
        k2=y1(j);
        k3=y1(j+1);
    else if (j==41)
            k1=y1(j-1);
            k2=y1(j);
            k3=y1(2);
        else
            k1=y1(j-1);
            k2=y1(j);
            k3=y1(j+1);
        end
    end
    
        k4= ((k3-k2)/dx1-(k2-k1)/dx1)/(2*dx1); % 2nd deriv
        r=((1+((k3-k1)/(2*dx1)).^2).^(3/2))/abs(k4);
        r1(j)=r;
   
end

c12=dx2/dx1;
c13=dx3/dx1;
c14=dx4/dx1;

for j=1:31 % both end points ignored
   if (j==1)
        k1=y2(30);
        k2=y2(j);
        k3=y2(j+1);
    else if (j==31)
            k1=y2(j-1);
            k2=y2(j);
            k3=y2(2);
        else
            k1=y2(j-1);
            k2=y2(j);
            k3=y2(j+1);
        end
    end
    
        k4= ((k3-k2)/dx2-(k2-k1)/dx2)/(2*dx2); % 2nd deriv
        r=((1+((k3-k1)/(2*dx2)).^2).^(3/2))/abs(k4);
        r2(j)=r;
   
end


for j=1:21 % both end points ignored
   if (j==1)
        k1=y3(20);
        k2=y3(j);
        k3=y3(j+1);
    else if (j==21)
            k1=y3(j-1);
            k2=y3(j);
            k3=y3(2);
        else
            k1=y3(j-1);
            k2=y3(j);
            k3=y3(j+1);
        end
    end
    
        k4= ((k3-k2)/dx3-(k2-k1)/dx3)/(2*dx3); % 2nd deriv
        r=((1+((k3-k1)/(2*dx3)).^2).^(3/2))/abs(k4);
        r3(j)=r;
   
end

for j=1:10001 % both end points ignored
   if (j==1)
        k1=y4(10000);
        k2=y4(j);
        k3=y4(j+1);
    else if (j==10001)
            k1=y4(j-1);
            k2=y4(j);
            k3=y4(2);
        else
            k1=y4(j-1);
            k2=y4(j);
            k3=y4(j+1);
        end
    end
    
        k4= ((k3-k2)/dx4-(k2-k1)/dx4)/(2*dx4); % 2nd deriv
        r=((1+((k3-k1)/(2*dx4)).^2).^(3/2))/abs(k4);
        dx41=dx4/100
        
        r4(j)=r;
   
end

figure(1); clf;
plot(x1,r1,'o-',x2,r2,'*-',x3,r3,'+-',x4,r4,'-')

figure(2); clf;
plot(x1,r1,'o-',x3,r3,'+-')

figure(3); clf;
plot(x4,r4,'+-')

Subject: curvature radius

From: Junghun

Date: 16 Nov, 2010 03:19:04

Message: 9 of 9

"Junghun " <jcho7@wisc.edu> wrote in message <ibqcff$563$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i9qctk$cat$1@fred.mathworks.com>...
> > "Junghun " <jcho7@wisc.edu> wrote in message <i9oes7$na4$1@fred.mathworks.com>...
> > > .........
> > > k1(1)=x1(j-1);
> > > k1(1)=x1(j);
> > > k1(1)=x1(j+1);
> > > .........
> > - - - - - - - - - - -
> > If the above is the code you actually used, you might have seen some strange results. Notice that for 1 < j < 21, k2(1) and k3(1) remain unchanged from the values they had at j = 1. I am sure this is not what you intended. Didn't you mean this:?
> >
> > k1(1)=x1(j-1);
> > k2(1)=x1(j);
> > k3(1)=x1(j+1);
> >
> > Just for the "record", for finding both a circle's center and radius the following requires fewer operations than the code I gave earlier. Again the three given points on the circle are (x1,y1), (x2,y2), and (x3,y3).
> >
> > x21 = x2-x1; y21 = y2-y1;
> > x31 = x3-x1; y31 = y3-y1;
> > h21 = x21^2+y21^2; h31 = x31^2+y31^2;
> > d = 2*(x21*y31-x31*y21);
> > a = x1+(h21*y31-h31*y21)/d;
> > b = y1-(h21*x31-h31*x21)/d;
> > r = sqrt(h21*h31*((x3-x2)^2+(y3-y2)^2))/abs(d);
> >
> > The circle's center is at (a,b) and its radius is r.
> >
> > Roger Stafford
>
> Thanks Roger,
> I am having same kind of problems that I asked before.
>
> Now, a curve is given. say, y=5.05*10.^-8*sin(10.^4*x)+0.4302*10.^-6
> I hope to get same curvature radius, r graph even if I change dx scale.
>
> I attach my code for that. (which is rough.)
> In the code below, I set 4 different x intervals, and got the r graph.
> The curvature radius is denpending on dx. (which I don't want to get)
>
> Even if I tried the vector calculus that you taught me,
> (by finding the radius of a circle which includes 3 given points, with vector calculus)
> r is still depending on dx.
>
> In the fixed region, I hope to get the same r graph even when dx is different.
> Could you help me out with this?
>
>
>
> =======================================
> clear all;
>
> x1 =(6*pi)/(40*10.^4):(6*pi)/(40*10.^4):(6*pi)/(40*10.^4)*41;
> x2 =(6*pi)/(30*10.^4):(6*pi)/(30*10.^4):(6*pi)/(30*10.^4)*31;
> x3 =(6*pi)/(20*10.^4):(6*pi)/(20*10.^4):(6*pi)/(20*10.^4)*21;
> x4 =(6*pi)/(10000*10.^4):(6*pi)/(10000*10.^4):(6*pi)/(10000*10.^4)*10001;
>
> y1 = 5.05*10.^-8*sin(10.^4*x1)+0.4302*10.^-6;
> y2 = 5.05*10.^-8*sin(10.^4*x2)+0.4302*10.^-6;
> y3 = 5.05*10.^-8*sin(10.^4*x3)+0.4302*10.^-6;
> y4 = 5.05*10.^-8*sin(10.^4*x4)+0.4302*10.^-6;
>
> dx1 =(6*pi)/(40*10.^4)
> dx2 =(6*pi)/(30*10.^4)
> dx3 =(6*pi)/(20*10.^4)
> dx4 =(6*pi)/(10000*10.^4)
>
> for j=1:41 % both end points ignored
> if (j==1)
> k1=y1(40);
> k2=y1(j);
> k3=y1(j+1);
> else if (j==41)
> k1=y1(j-1);
> k2=y1(j);
> k3=y1(2);
> else
> k1=y1(j-1);
> k2=y1(j);
> k3=y1(j+1);
> end
> end
>
> k4= ((k3-k2)/dx1-(k2-k1)/dx1)/(2*dx1); % 2nd deriv
> r=((1+((k3-k1)/(2*dx1)).^2).^(3/2))/abs(k4);
> r1(j)=r;
>
> end
>
> c12=dx2/dx1;
> c13=dx3/dx1;
> c14=dx4/dx1;
>
> for j=1:31 % both end points ignored
> if (j==1)
> k1=y2(30);
> k2=y2(j);
> k3=y2(j+1);
> else if (j==31)
> k1=y2(j-1);
> k2=y2(j);
> k3=y2(2);
> else
> k1=y2(j-1);
> k2=y2(j);
> k3=y2(j+1);
> end
> end
>
> k4= ((k3-k2)/dx2-(k2-k1)/dx2)/(2*dx2); % 2nd deriv
> r=((1+((k3-k1)/(2*dx2)).^2).^(3/2))/abs(k4);
> r2(j)=r;
>
> end
>
>
> for j=1:21 % both end points ignored
> if (j==1)
> k1=y3(20);
> k2=y3(j);
> k3=y3(j+1);
> else if (j==21)
> k1=y3(j-1);
> k2=y3(j);
> k3=y3(2);
> else
> k1=y3(j-1);
> k2=y3(j);
> k3=y3(j+1);
> end
> end
>
> k4= ((k3-k2)/dx3-(k2-k1)/dx3)/(2*dx3); % 2nd deriv
> r=((1+((k3-k1)/(2*dx3)).^2).^(3/2))/abs(k4);
> r3(j)=r;
>
> end
>
> for j=1:10001 % both end points ignored
> if (j==1)
> k1=y4(10000);
> k2=y4(j);
> k3=y4(j+1);
> else if (j==10001)
> k1=y4(j-1);
> k2=y4(j);
> k3=y4(2);
> else
> k1=y4(j-1);
> k2=y4(j);
> k3=y4(j+1);
> end
> end
>
> k4= ((k3-k2)/dx4-(k2-k1)/dx4)/(2*dx4); % 2nd deriv
> r=((1+((k3-k1)/(2*dx4)).^2).^(3/2))/abs(k4);
> dx41=dx4/100
>
> r4(j)=r;
>
> end
>
> figure(1); clf;
> plot(x1,r1,'o-',x2,r2,'*-',x3,r3,'+-',x4,r4,'-')
>
> figure(2); clf;
> plot(x1,r1,'o-',x3,r3,'+-')
>
> figure(3); clf;
> plot(x4,r4,'+-')


Sorry, I figured it out. I've done some stupid things.

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