"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(j1);
> > > 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(j1);
> > 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 = x2x1; y21 = y2y1;
> > x31 = x3x1; y31 = y3y1;
> > h21 = x21^2+y21^2; h31 = x31^2+y31^2;
> > d = 2*(x21*y31x31*y21);
> > a = x1+(h21*y31h31*y21)/d;
> > b = y1(h21*x31h31*x21)/d;
> > r = sqrt(h21*h31*((x3x2)^2+(y3y2)^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(j1);
> k2=y1(j);
> k3=y1(2);
> else
> k1=y1(j1);
> k2=y1(j);
> k3=y1(j+1);
> end
> end
>
> k4= ((k3k2)/dx1(k2k1)/dx1)/(2*dx1); % 2nd deriv
> r=((1+((k3k1)/(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(j1);
> k2=y2(j);
> k3=y2(2);
> else
> k1=y2(j1);
> k2=y2(j);
> k3=y2(j+1);
> end
> end
>
> k4= ((k3k2)/dx2(k2k1)/dx2)/(2*dx2); % 2nd deriv
> r=((1+((k3k1)/(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(j1);
> k2=y3(j);
> k3=y3(2);
> else
> k1=y3(j1);
> k2=y3(j);
> k3=y3(j+1);
> end
> end
>
> k4= ((k3k2)/dx3(k2k1)/dx3)/(2*dx3); % 2nd deriv
> r=((1+((k3k1)/(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(j1);
> k2=y4(j);
> k3=y4(2);
> else
> k1=y4(j1);
> k2=y4(j);
> k3=y4(j+1);
> end
> end
>
> k4= ((k3k2)/dx4(k2k1)/dx4)/(2*dx4); % 2nd deriv
> r=((1+((k3k1)/(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.
