求助：最优曲线拟合（对数螺旋线）

lhq

已知曲线由以下方程组确定：

(x+2.5)^2+(y+4.8)^2=25^2

x=-12

x=12

y>0

根据以上条件，可知曲线是在（-12<x<12）且（y>0）之间的一段圆弧线！

现需要寻找一段对数螺旋线，要求其和已知的曲线最接近，也就是最优吧！

想用ＭＡＴＬＡＢ软件做曲线拟合，但苦于能力有限，故求助各位高手！！！

对数螺旋线方程：r=r0*exp^(m*k)       (极坐标)

注：参数r0  是初始极径

       参数m  是极角

       参数k   是渐开线系数

希望高手们答复，最好有编好的程序！！！感激不尽！

mhkmars

1. x=linspace(-12,12,100);
   
2. y=sqrt(25^2-(x+2.5).^2)-4.8;
   
3. plot(x,y,'b*')
   
4. axis equal
   
5. [th,r]=cart2pol(x,y);
   
6. [para,res]=lsqcurvefit(@myfun,[10,2],th,r);
   
7. cal_r=para(1)*exp(para(2)*th);
   
8. [cal_x,cal_y]=pol2cart(th,cal_r);
   
9. hold on
   
10. plot(cal_x,cal_y,'r^')

复制代码

lhq

不局限于ＭＡＴＬＡＢ软件，也可用其它软件！

lhq

mhkmars 发表于 2012-4-26 16:40

首先，非常感谢你的帮助！还有几个问题想向你请教一下：
1.拟合出来的对数螺旋线，左侧端点可以比原来的圆弧线上的左侧端点高一点吗？如果能，那是最好了！希望你帮我再考虑一下！
2.小弟刚学MATLAB不久，你发的程序，看不太明白！希望你能给我发点这方面的学习资料，小弟邮箱：lihaiquanemail@163.com，万分感激！！！

lhq

mhkmars 发表于 2012-4-26 16:40

3.能不能把拟合出来的对数螺旋线的几个参数也求出来呢？就是初始极径、极角和渐开线系数！

mhkmars

lhq 发表于 2012-4-27 16:35
首先，非常感谢你的帮助！还有几个问题想向你请教一下：
1.拟合出来的对数螺旋线，左侧端点可以比原来的 ...

1. 
   
2. >> para
   
3. 
   
4. para =
   
5. 
   
6.    17.5453    0.0935

复制代码

1.为什么要高一点，你这么做的目的是什么？matlab做拟合都是基于一定的收敛准则的，在你这些参数的设定下求的的最优值就是这个，所以没有改的必要。
2.matlab的资料网上一搜索一大把，在simwe上面也可以找到你想要的资料，另外，我认为最好的matlab的学习资料就是matlab 的help文档。你可以从这方面开始学。多写写代码，多上simwe上泡泡，版上还是有很多人很乐意帮助你的。
3.看来确实是学习不久(⊙o⊙)…我给你的代码里面的para就是你要求解的初始极径和你的系数
还是贴出来吧。。。

lhq

mhkmars 发表于 2012-4-27 18:34
1.为什么要高一点，你这么做的目的是什么？matlab做拟合都是基于一定的收敛准则的，在你这些参数的设定下 ...

非常感谢！

lhq

mhkmars 发表于 2012-4-26 09:55

你好，我用ＭＡＴＬＡＢ运行你发的程序，到第6步时，就提示有错误啊！你看看是不是发程序时有什么遗漏或错误？请你帮忙看看！谢了１

mhkmars

lhq 发表于 2012-4-28 10:46
你好，我用ＭＡＴＬＡＢ运行你发的程序，到第6步时，就提示有错误啊！你看看是不是发程序时有什么遗漏或 ...

1. function out=myfun(x,xdata)
   
2. out=x(1)*exp(x(2)*xdata);

复制代码

不好意思，忘记把myfun这个函数发给你了

lhq

mhkmars 发表于 2012-4-28 15:43
不好意思，忘记把myfun这个函数发给你了

麻烦你把整个程序的完整步骤再给我重发一下吧！要是能附带些步骤的解析就更感谢了，实在是刚学ＭＡＴＬＡＢ,有点儿吃力呀！所以还请你再帮帮忙啦！谢啦！

mhkmars

lhq 发表于 2012-4-28 20:01
麻烦你把整个程序的完整步骤再给我重发一下吧！要是能附带些步骤的解析就更感谢了，实在是刚学ＭＡＴＬＡ ...

myfun函数为：

1. function out=myfun(x,xdata)
   
2. %这个函数就是你要进行拟合的函数的形式
   
3. %复制代码后，保存为一个名字和函数名myfun同名的文件名
   
4. out=x(1)*exp(x(2)*xdata);

复制代码

1. x=linspace(-12,12,100);%从-12到12产生100个数
   
2. y=sqrt(25^2-(x+2.5).^2)-4.8;%计算对应于x的相应的y
   
3. plot(x,y,'b*')%绘制x-y的函数图
   
4. axis equal
   
5. [th,r]=cart2pol(x,y);%将直角坐标转换为极坐标
   
6. [para,res]=lsqcurvefit(@myfun,[10,2],th,r);%开始曲线拟合，para是你需要的参数，res是残差
   
7. cal_r=para(1)*exp(para(2)*th);%利用拟合后的参数重新计算幅角和半径
   
8. [cal_x,cal_y]=pol2cart(th,cal_r);%极坐标形式变换为直角坐标
   
9. hold on
   
10. plot(cal_x,cal_y,'r^')%绘制拟合后的图线

复制代码

lhq

mhkmars 发表于 2012-4-28 20:09
myfun函数为：

谢谢啊！我自己研究研究，呵呵！

lhq

mhkmars 发表于 2012-4-26 09:55

你好，向你请教几个问题！
这是你给我的ＭＡＴＬＡＢ程序：
function out=myfun(x,xdata)
out=x(1)*exp(x(2)*xdata);

x=linspace(-12,12,100);
y=sqrt(25^2-(x+2.5).^2)-4.8;
plot(x,y,'b*')
axis equal
[th,r]=cart2pol(x,y);
[para,res]=lsqcurvefit(@myfun,[10,2],th,r);
cal_r=para(1)*exp(para(2)*th);
[cal_x,cal_y]=pol2cart(th,cal_r);
hold on
plot(cal_x,cal_y,'r^')

1.程序中的[10，2]是什么意思？
2.你的程序我已经理解了，你是用已知曲线的直坐标点转化为相应的极坐标点，从而实现对数螺旋线的拟合！这就相当于已经限定了对数螺旋线的极角范围就是从（52度到123度之间吧）！从ＭＡＴＬＡＢ拟合出来的图形看，在限定极角的条件下，拟合出来的曲线是这样，但如果不把对数螺旋线的极角限定在这个范围，而是从两条曲线的拟合度来考虑，相信肯定存在一条更接近的对数螺旋线，对吧？这么来理解，从你拟合出来的曲线来看，如果在对数螺旋线的左侧再取一段延长线的话，那它肯定比（52度到123度）这段更接近原来的圆弧曲线！对吧？这一点是肯定的！不知，我的想法你了解没有？
    希望能和你交流学习！

lhq

我在想可不可以按这种思路分析，把对数螺旋线的极角设定成一个变化值，即（x+52,x+123），但不知这种思路可否用ＭＡＴＬＡＢ实现？

mhkmars

1.[10 2]是指优化的初始值。
2.你提到的52到123我是从你的给的x的范围是从-12到12的来的。因为你已经给我了一个待拟合的模型，如果你想进一步提高拟合的效果，可以考虑修改你的模型。

lhq

mhkmars 发表于 2012-5-2 11:38
1.[10 2]是指优化的初始值。
2.你提到的52到123我是从你的给的x的范围是从-12到12的来的。因为你已经给我了 ...

1.　我知道[10，2]是初始值，但我不知道为什么这么取？为什么要取10，为什么要取2？
2.　我明白你说的，是按给定的直坐标转化的，也就是说限制了对数螺旋线的极角范围！但现在为了找到更好的拟合曲线，不限制它的极角范围，把它的极角范围看成是（x+52,x+123），可以把x也看成是一个参数，这样的意思！你从这个思路出发，看看能不能把程序改改，以实现这种效果？
3.　para输出来的极角，也就是th是弧度制，能不能在程序中就改成角度制，让它输出成角度制，因为这个要影响到para输出的另一个参数！拿实际数据来说，程序输出的是para [17.5452，0.0935]，th1=2.1506，但经计算验证，如果把th1转化成2.1506*180/pi，那么对应的para　应该是　[17.5452，0.001632054]　才对！！！

希望你帮我考虑考虑以上3个问题，感激不尽！

mhkmars

lhq 发表于 2012-5-2 16:34
1.　我知道[10，2]是初始值，但我不知道为什么这么取？为什么要取10，为什么要取2？
2.　我明白你说的， ...

1.初始值是我自己试着给的，最后反正能够收敛就行，至于为什么，你可以看看相关的数值分析的书
2.没怎么看懂你的。。。但是你可以试着修改一下你的拟合模型，其实效果更好。我将你的模型后面加上一个常数项，变为r=r0*exp（k*theta）+D；拟合的误差要小很多；具体如下
修改前的拟合误差：

1. res =
   
2. 
   
3.     0.1387

复制代码

修改后的拟合误差：

1. res =
   
2. 
   
3.     0.0013

复制代码

代码修改了一下
模型函数为：

1. function out=myfun(x,xdata)
   
2. %这个函数就是你要进行拟合的函数的形式
   
3. %复制代码后，保存为一个名字和函数名myfun同名的文件名
   
4. out=x(1)*exp(x(2)*xdata)+x(3);

复制代码

运行下面的程序段

1. x=linspace(-5,5,100);%从-12到12产生100个数
   
2. y=sqrt(25^2-(x+2.5).^2)-4.8;%计算对应于x的相应的y
   
3. plot(x,y,'b*')%绘制x-y的函数图
   
4. axis equal
   
5. [th,r]=cart2pol(x,y);%将直角坐标转换为极坐标
   
6. [para,res]=lsqcurvefit(@myfun,[10,2,0],th,r);%开始曲线拟合，para是你需要的参数，res是残差
   
7. cal_r=para(1)*exp(para(2)*th)+para(3);%利用拟合后的参数重新计算幅角和半径
   
8. [cal_x,cal_y]=pol2cart(th,cal_r);%极坐标形式变换为直角坐标
   
9. hold on
   
10. plot(cal_x,cal_y,'r^')%绘制拟合后的图线

复制代码

拟合的图像


3.没看懂th1是个什么。。。
以上！

lhq

mhkmars 发表于 2012-5-13 19:08
1.初始值是我自己试着给的，最后反正能够收敛就行，至于为什么，你可以看看相关的数值分析的书
2.没怎么 ...

你好，对于你的交流学习，本人表示非常感谢！
我按照你的建议试着研究了一下，觉得拟合效果可以，但感觉不是最理想！
所以，我又考虑能不能从相对坐标系的角度出发，去研究这个问题？
我结合附件中的图，说明一下我的思想，我觉得从数学理论的角度肯定是可行，但在MATLAB中怎么转化，我还需要你的帮忙啊！

     假设图中的曲线MDH就是我们所要拟合的曲线，结合给出的方程我们知道它是一段以O1为圆心的圆弧，X范围介于（-12，12）之间！这段曲线其实是可求的，包括其上的各点！
我现在以最右端的端点M点为原点建立一个相对极角坐标系，即XMY；在这个极坐标系中，求出曲线MDH上的一系列点值，即A、B、C、D、E……H点相对于原点M的极角和极径，利用
这些求得的点在相对坐标系中寻找一段和其最相似的对数螺旋线，使寻找到的那段对数螺旋线和已求的一系列点的拟合度最高！最后，再把拟合出来的对数螺旋线转化到原来的笛卡尔
坐标系中！基本思想就是这样！
      朋友，希望你有时间帮我看看，谢谢了！
    如果方便的话，留个电话，我们可以电话沟通！（本人：15242241758）