1stopt如何求解超定方程组

szldh2005

设有方程组：
   x'=a1*x+b1*y+c1
   y'=a2*x+b2*y+c2
今有n组（x,y）及（x',y'）数据，请教，如何用1stopt求解参数a1,b1,c1,a2,b2,c2?
谢谢！

shamohu

等同于微分方程拟合，1stOpt 4.0版开始有此功能，楼主的问题应该可以很容易解决。方便的话把数据贴上来。

szldh2005

Title "超定方程组的系数求解";
//超定方程组的系数求解，是用共享模式求解吗？pro1.5
Parameters ix,iy,k,a ;
Variable xs,ys,x3,y3 ;
//sharedModel;   //可以用共享模式吗？
Function x3=ix+k*cos(a)*xs+k*sin(a)*ys;
         y3=iy-k*sin(a)*xs+k*cos(a)*ys;
data;

shamohu

因为参数a、k两个公式中都有，因此必须用“SharedModel”。
结果几乎完全重合，楼主是自己赋值造的一组数据吧！

均方差(RMSE): 0.000559220214480197
残差平方和(SSE): 4.19054512699592E-5
相关系数(R): 1
相关系数之平方(R^2): 1

参数                  最佳估算
--------------------        -------------
ix        2472721.34671741
iy        391032.172340117
k        0.999991068736129
a        0.0170804450845957

wanglu

应该有好多解吧？
以下是多次运行Forcal代码的结果：

1. !using["fcopt","math","sys"];
   
2. init(::Array,max)=
   
3. {
   
4.     max=67,
   
5.     Array=arrayinitns{max,4 :
   
6.         "21931.951  117563.229  2496657.820  508202.617
   
7. 14819.722  99053.090  2489230.549  489816.816
   
8. 49990.336  87182.164  2524192.968  477347.032
   
9. 41838.752  140307.803  2516950.004  530603.673
   
10. 20258.393  157719.614  2495670.370  548381.370
    
11. 26548.263  106139.151  2501078.300  496701.464
    
12. 31825.976  104925.249  2506334.463  495397.608
    
13. 34434.524  99218.909  2508845.145  489647.599
    
14. 25467.204  92605.317  2499766.259  483188.188
    
15. 27506.217  96237.631  2501866.993  486785.114
    
16. 34602.641  90096.720  2508857.436  480523.951
    
17. 35575.694  94853.603  2509911.583  485263.479
    
18. 37645.160  103099.408  2512121.563  493472.662
    
19. 42774.324  112412.213  2517408.989  502696.422
    
20. 43210.824  101659.490  2517661.770  491937.908
    
21. 38334.616  107930.793  2512893.428  498291.523
    
22. 41648.750  110013.427  2516242.619  500317.231
    
23. 37351.833  114320.444  2512019.929  504696.970
    
24. 40794.893  96636.575  2515160.425  486957.034
    
25. 49291.723  91964.693  2523576.146  482140.754
    
26. 45833.955  90167.433  2520088.218  480402.828
    
27. 51036.764  98679.130  2525435.597  488824.347
    
28. 45881.631  93853.881  2520198.849  484087.891
    
29. 46285.362  99141.362  2520692.825  489367.657
    
30. 48495.045  104124.991  2522987.283  494312.776
    
31. 45695.851  85121.477  2519863.953  475360.012
    
32. 40277.662  89711.063  2514524.991  480041.427
    
33. 30775.459  120254.569  2505545.925  510742.498
    
34. 44993.803  113077.243  2519639.483  503323.441
    
35. 36779.076  124432.038  2511619.961  514816.782
    
36. 38004.330  129249.125  2512927.299  519612.196
    
37. 32964.556  129215.303  2507887.728  519664.456
    
38. 49187.733  132799.223  2524169.605  522970.739
    
39. 44064.522  135740.635  2519097.424  525999.197
    
40. 40758.132  131890.827  2515725.794  522206.456
    
41. 34494.553  134958.674  2509515.582  525380.806
    
42. 23122.600  107379.121  2497674.346  497999.750
    
43. 15908.879  114823.693  2490588.890  505566.376
    
44. 25487.773  111230.576  2500104.933  501810.213
    
45. 20092.921  118566.578  2494836.212  509237.221
    
46. 33843.010  153045.304  2509173.049  543475.765
    
47. 37632.144  154282.492  2512982.726  544648.045
    
48. 27803.508  164725.056  2503333.966  555256.861
    
49. 28149.018  158645.975  2503575.595  549172.820
    
50. 20917.019  163613.134  2496429.552  554262.728
    
51. 16044.041  169303.661  2491654.520  560035.602
    
52. 12870.303  166434.768  2488432.274  557221.359
    
53. 36296.931  145677.097  2511500.745  536066.787
    
54. 28779.827  153445.273  2504117.482  543962.149
    
55. 27098.110  153300.973  2502433.560  543846.594
    
56. 24261.566  159486.729  2499703.105  550079.839
    
57. 30516.857  144860.079  2505707.612  535348.616
    
58. 18275.734  160873.067  2493741.878  551568.196
    
59. 9857.462  161992.356  2485344.025  552831.093
    
60. 13399.860  161233.489  2488872.914  552011.841
    
61. 34456.507  140344.377  2509569.526  530766.325
    
62. 44540.518  144398.723  2519721.222  534647.814
    
63. 49314.448  144872.132  2524502.498  535039.614
    
64. 39505.030  149609.090  2514775.504  539943.379
    
65. 47828.309  139148.110  2522918.827  529341.861
    
66. 27183.099  126133.035  2502054.522  516681.409
    
67. 19595.538  123807.305  2494428.414  514485.631
    
68. 23394.983  131399.060  2498356.934  522011.318
    
69. 27555.024  144141.015  2502733.956  534680.249
    
70. 24306.484  139503.599  2499406.715  530099.036
    
71. 31487.943  136831.141  2506541.417  527304.336
    
72. 25987.014  117614.073  2500713.124  508184.195"
    
73.     }.free()
    
74. };
    
75. f(ix, iy, k, a :i,s,xs,ys,x3,y3:Array,max)=
    
76. {
    
77.     s=0,i=0,(i<max).while{
    
78.         Array.GA[i*4, &xs,&ys,&x3,&y3],
    
79.         s=s+[ ix+k*cos(a)*xs+k*sin(a)*ys-x3]^2+[ iy-k*sin(a)*xs+k*cos(a)*ys-y3]^2,
    
80.         i++
    
81.     },
    
82.     sqrt[s/max]
    
83. };
    
84. Opt[HFor("f")];

复制代码

结果：
2472721.346717395         391032.1723404035         0.9999910687332905        -12.54929016927396        7.908568261468204e-004
2472721.346717022         391032.172339913          0.9999910687382372        -106.7970697769665        7.908567985651832e-004
2472721.346717131         391032.1723400778         0.9999910687370118        1.7080445085375e-002      7.908567901771496e-004
2472721.346717442         391032.1723402201         0.9999910687353089        -6.266104862095189        7.908568057345164e-004

wanglu

还有：
2472721.346717457         391032.1723403272         0.9999910687341977        14677.5379580166          5.592202426691267e-004
2472721.345640144         391032.1726446233         0.9999910687338469        -741010651.0432801        5.826017927098407e-004

确实有很多解。

szldh2005

shamohu版主：
我给出的pro1.5版本的代码程序是否正确？因为在pro1.5上运行2h都没有得到结果。数据是实际资料。
版主运行的版本号，运行时间大概是多少？
版主用1stopt解算的结果及wanglu网友用Forca求解的第3组解满足工程实际。
谢谢网友给予的大力帮助。祝网友新年万事顺意!

shamohu

将“//sharedModel;”改为“sharedModel;”就行了，其它都正确。

用的4.0版，计算时间也就几秒。

szldh2005

shamohu版主：
再请教一下，这个例子，pro1.5版本是不是不能求解，最低版本是要多少？谢谢！