微分方程参数优化（拟合）问题

wanglu

在Matlabsky论坛和Matlab中文论坛有两个微分方程参数拟合问题，我用OpenFC（Forcal程序）求解了，但没有得到验证，希望这里的高手能帮忙验证一下，同时讨论一下此类问题的不同解法。

OpenFC下载： http://xiazai.zol.com.cn/detail/27/262791.shtml

首先说明一下，对于参数拟合，我用的是徐士良算法库中的求n维极值的单形调优法函数XSLSF::jsim，该函数需要试算求解，有时不能获得最优解。对于这两个例子，第一个自我感觉良好，似乎得到了最优解；但第二个运算量太大，十几分钟才获得一组解，难以尝试各种参数，故仅获得了迭代100次的局部解。

第一个例子地址：http://www.matlabsky.cn/thread-3511-1-3.html

据说本例可以用Matlab的dsolve函数得出y和z的表达式，因而验证较容易，麻烦Matlab高手能给出解析解。

数学模型：

1. 
   
2.       dy/dt=k*y*z+0.095*b*z
   
3.       dz/dt=-b*z-0.222*z
   

复制代码

实验数据（ti; yi）：
ti        yi
0.01    2.329101563
0.68    3.851292783
1.1     4.50093936
1.63    6.749172247
2.07    9.112062872
2.67    9.691657366
3.09    11.16928013
3.64    10.91451823
4.65    16.44279204
5.1     18.29615885
5.58    21.63989025
6.11    25.78611421
6.63    26.34282633
7.06    26.50581287
7.62    27.63951823
8.66    35.02757626
9.04    35.5562035
9.63    36.10393415

    已知z在t=0.01时的初值为5000，求k和b的最优值？

    算法分析：用求n维极值的单形调优法求最优参数值k和b。用连分式法对微分方程组积分一步函数pbs1求理论值，与实验值做比较，理论值与实验值差的平方和最小为最优。

    Forcal代码：

1. 
   
2. !using["XSLSF"];                //使用命名空间XSLSF
   
3. //函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。
   
4. //t为自变量，y和z为函数值，dy和dz为右端函数值（即微分方程的值）。
   
5. //该函数被调用时将用到参数k和b，这2个参数正好是拟合参数，用模块变量传递这2个参数。
   
6. f(t,y,z,dy,dz::k,b)=
   
7. {
   
8.   dy=k*y*z+0.095*b*z,
   
9.   dz=-b*z-0.222*z
   
10. };
    
11. //用连分式法对微分方程组进行积分，获得理论值。
    
12. //t1,t2为积分的起点和终点。
    
13. //h,s为自动变量。
    
14. //模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。
    
15. 积分(t1,t2:h,s:hf,Array,step,eps)=
    
16. {
    
17.     s=ceil[(t2-t1)/step],        //计算积分步数
    
18.     h=(t2-t1)/s,                 //重新计算积分步长
    
19.     {   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
    
20.         t1=t1+h                  //积分步长增加
    
21.     }.until[abs(t1-t2)<h/2]      //连续积分至t2
    
22. };
    
23. //目标函数定义，自变量_k,_b为需要优化的参数，需要将这些参数传递给对应的模块变量k,b。
    
24. //模块变量：Array为工作数组；数组tArray存放实验数据ti；数组yArray存放实验数据yi；max为数组tArray和yArray的长度；k,b为优化变量。
    
25. 优化(_k,_b:i,s,y,z,yy,t1,t2:Array,tArray,yArray,max,k,b)=
    
26. {
    
27.   k=_k,b=_b,                         //传递优化变量，函数f中要用到k,b
    
28.   Array.setra(0,2.329101563,5000),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
29.   i=1,s=0,t1=0.01,
    
30.   (i<max).while{
    
31.     tArray.getra(i,&t2), yArray.getra(i,&yy),
    
32.     积分(&t1,t2),                    //从t1积分到t2
    
33.     Array.getra(0,&y,z),             //从数组Array获得t2时的积分值
    
34.     s=s+[y-yy]^2,                    //累加理论值与实验值差的平方
    
35.     i++
    
36.   },
    
37.   s
    
38. };
    
39. 验证(_k,_b:i,y,z,yy,t1,t2:Array,tArray,yArray,max,k,b)=
    
40. {
    
41.   k=_k,b=_b,
    
42.   printff{"\r\n               t            y实验值            y模拟值         z模拟值\r\n"},
    
43.   Array.setra(0,2.329101563,5000),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
44.   i=1,t1=0.01,
    
45.   (i<max).while{
    
46.     tArray.getra(i,&t2), yArray.getra(i,&yy),
    
47.     积分(&t1,t2),                    //从t1积分到t2
    
48.     Array.getra(0,&y,&z),            //从数组Array获得t2时的积分值
    
49.     printff{"{1,r,18.8}{2,r,18.8}{3,r,18.8}{4,r,18.8}\r\n",t1,yy,y,z},
    
50.     i++
    
51.   }
    
52. };
    
53. main(:x,xx,g,d,u,v,k,i,t1,t2,t3,_eps:Array,tArray,yArray,max,hf,step,eps)=
    
54. {
    
55.   hf=HFor("f"),                                    //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
    
56.   step=0.01,eps=1e-6,                              //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
    
57.   Array=new[rtoi(real_s),rtoi(2*15)],              //申请工作数组
    
58.   max=18,                                          //实验数据组数
    
59.   tArray=new{rtoi(real_s),rtoi(max),rtoi(EndType), //存放实验数据ti
    
60.     0.01,
    
61.     0.68,
    
62.     1.1,
    
63.     1.63,
    
64.     2.07,
    
65.     2.67,
    
66.     3.09,
    
67.     3.64,
    
68.     4.65,
    
69.     5.1,
    
70.     5.58,
    
71.     6.11,
    
72.     6.63,
    
73.     7.06,
    
74.     7.62,
    
75.     8.66,
    
76.     9.04,
    
77.     9.63
    
78.     },
    
79.   yArray=new{rtoi(real_s),rtoi(max),rtoi(EndType), //存放实验数据yi
    
80.     2.329101563,
    
81.     3.851292783,
    
82.     4.50093936,
    
83.     6.749172247,
    
84.     9.112062872,
    
85.     9.691657366,
    
86.     11.16928013,
    
87.     10.91451823,
    
88.     16.44279204,
    
89.     18.29615885,
    
90.     21.63989025,
    
91.     25.78611421,
    
92.     26.34282633,
    
93.     26.50581287,
    
94.     27.63951823,
    
95.     35.02757626,
    
96.     35.5562035,
    
97.     36.10393415
    
98.     },
    
99.   x=new[rtoi(real_s),rtoi(3)],                     //申请工作数组
    
100.   xx=new[rtoi(real_s),rtoi(2),rtoi(3)],            //申请工作数组
     
101.   g=new[rtoi(real_s),rtoi(3)],                     //申请工作数组
     
102.   _eps=1e-100, d=550,u=1.6,v=0.4,k=500,            //变换d、u、v进一步求解，k为允许的最大迭代次数
     
103.   i=jsim[HFor("优化"),d,u,v,x,_eps,k,xx,g],        //求n维极值的单形调优法
     
104.   x.getra(0,&t1,&t2,&t3),                          //获得最优参数值及目标终值
     
105.   printff{"\r\n实际迭代次数={1,i}, k={2,r}, b={3,r}, 目标函数终值={4,r}\r\n",i,t1,t2,t3},
     
106.   验证(t1,t2),                                     //验证计算结果
     
107.   delete[Array],delete[tArray],delete[yArray],delete[x],delete[xx],delete[g]
     
108. };
     

复制代码

结果：

1. 
   
2. 实际迭代次数=104, k=1.4083232846396376e-004, b=-1.4666335917093112e-004, 目标函数终值=24.291246634355627
   
3. 
   
4.                t            y实验值            y模拟值         z模拟值
   
5.               0.68         3.8512928         3.5561928         4309.3892
   
6.                1.1         4.5009394         4.5088671          3925.987
   
7.               1.63         6.7491722         5.9132744         3490.4669
   
8.               2.07         9.1120629         7.2438212         3165.8451
   
9.               2.67         9.6916574          9.277481         2771.2765
   
10.               3.09          11.16928         10.832584         2524.7187
    
11.               3.64         10.914518         13.002227         2234.7074
    
12.               4.65         16.442792         17.253317         1786.1071
    
13.                5.1         18.296159         19.204517         1616.4048
    
14.               5.58          21.63989         21.291108         1453.1227
    
15.               6.11         25.786114         23.574516          1291.924
    
16.               6.63         26.342826         25.768822         1151.1585
    
17.               7.06         26.505813          27.53363         1046.4169
    
18.               7.62         27.639518         29.747528         924.16377
    
19.               8.66         35.027576         33.556979          733.7452
    
20.               9.04         35.556204         34.840955         674.42303
    
21.               9.63         36.103934         36.714588          591.6789
    

复制代码

wanglu

第二个问题，将允许的最大迭代次数k＝500，运行了约40分钟，终于获得一组最优解：

1. 
   
2. 实际迭代次数=213, k1=9.5111080627706936, k2=18.593155843279298, k4=-0.3629342128423555, k11=9.2785258647078628, k16=42.278370674583023, k17=23.125269056444612, 目标函数终值=1.0494952055839093e-009
   

复制代码

另外，第一个问题有解析解的，Matlab、maple或mmtc高手们能否给出解析解呢？

盼望大家帮忙验证一下我的结果，给出目标函数终值，最好能给出代码，让大家都能参考。

欢迎大家讨论一下此类问题的解法，无论是用Matlab、maple、mmtc或1stOpt，还是用其他软件均可。

wanglu

第二个例子：http://www.ilovematlab.cn/thread-57183-1-4.html

现有如下6个微分式，k1, k2, k4, k11, k16, k17为待定的常数，与t无关，目标是通过[A]、[B]、[C]、[D]、[E]、[F]关于t的函数解析式，并用实验数据（[E]关于时间t的一组数据），拟合求出这些待定量k的值，实验数据见excel文件。
如果求不出函数解析解的话，也可用其他方法把待定的常数k值求出来也行。

已知t=0时，[A]=[C]=[F]=0, [B]=1.25*10-3, [D]= 1.0*10-4, [E]= 2.0*10-4

d[A]/dt=56*k1*[D]-k2*[A]*[B]
d[B]/dt=-k2*[A]*[B]
d[C]/dt=k2*[A]*[B]-k4*[C]*[D]-k16*[C]*[E]-56*k17*[C]*[F]
d[D]/dt=0.002*k11*[A]- 56*k1*[D]
d[E]/dt=-k16*[C]*[E]
d[F]/dt= k16*[C]*[E]-56*k17*[C]*[F]

t/s	E浓度/(mg/L)
0	0.0002
1	0.000199566
2	0.000199814
3	0.000199068
4	0.000199006
5	0.000198448
6	0.000198634
7	0.000197704
8	0.000197642
9	0.000197454
10	0.000197082
11	0.00019671
12	0.0001964
13	0.000196276
14	0.000195716
15	0.000195468
16	0.000195034
17	0.000194724
18	0.000194166
19	0.000193978
20	0.00019342
21	0.00019311
22	0.000192738
23	0.000192428
24	0.000192054
25	0.000191558
26	0.000191124
27	0.000191
28	0.00019044
29	0.000189944
30	0.000189882
31	0.000189076
32	0.000188702
33	0.000188268
34	0.000187958
35	0.000187524
36	0.000187088
37	0.000186654
38	0.000186096
39	0.000185724
40	0.000185164
41	0.000184668
42	0.000184358
43	0.00018386
44	0.00018355
45	0.000183116
46	0.000182496
47	0.000181998
48	0.000181502
49	0.000181068
50	0.000180696
51	0.000180322
52	0.000179578
53	0.000179268
54	0.000178584
55	0.000178212
56	0.000177654
57	0.000177094
58	0.000176722
59	0.000176226
60	0.000175606
61	0.00017517
62	0.00017455
63	0.00017424
64	0.00017368
65	0.000172998
66	0.000172564
67	0.000171942
68	0.000171508
69	0.000171074
70	0.000170392
71	0.000169956
72	0.000169336
73	0.000168778
74	0.00016828
75	0.00016766
76	0.000167288
77	0.000166728
78	0.000166108
79	0.00016555
80	0.000164928
81	0.000164494
82	0.000163998
83	0.000163376
84	0.00016288
85	0.00016226
86	0.000161824
87	0.000161204
88	0.000160584
89	0.000160024
90	0.000159528
91	0.000159032
92	0.00015841
93	0.000157852
94	0.000157232
95	0.000156672
96	0.000156176
97	0.000155556
98	0.000155122
99	0.000154438
100	0.00015388
101	0.00015332
102	0.000152762
103	0.000152204
104	0.000151582
105	0.000151086
106	0.000150528
107	0.000149906
108	0.00014941
109	0.00014879
110	0.00014823
111	0.00014761
112	0.000147114
113	0.00014643
114	0.000145996
115	0.000145376
116	0.000144816
117	0.000144258
118	0.0001437
119	0.000143078
120	0.000142458
121	0.000141962
122	0.000141402
123	0.000140844
124	0.000140224
125	0.000139664
126	0.000139106
127	0.000138548
128	0.000137988
129	0.00013743
130	0.00013681
131	0.000136312
132	0.000135692
133	0.000135134
134	0.000134512
135	0.000134016
136	0.000133396
137	0.000132836
138	0.000132278
139	0.00013172
140	0.00013116
141	0.000130602
142	0.000130106
143	0.000129484
144	0.000128864
145	0.000128368
146	0.000127808
147	0.00012725
148	0.000126568
149	0.000126132
150	0.000125574
151	0.000124954
152	0.000124456
153	0.000123836
154	0.000123278
155	0.00012278
156	0.000122222
157	0.000121602
158	0.000121104
159	0.000120484
160	0.000119988
161	0.000119428
162	0.00011887
163	0.000118312
164	0.000117752
165	0.000117194
166	0.000116698
167	0.00011614
168	0.00011558
169	0.000115084
170	0.000114526
171	0.000113966
172	0.000113408
173	0.00011285
174	0.00011229
175	0.000111794
176	0.000111236
177	0.000110676
178	0.000110118
179	0.000109622
180	0.000109124
181	0.000108566
182	0.000108008
183	0.000107448
184	0.00010689
185	0.000106394
186	0.000105834
187	0.000105338
188	0.00010478
189	0.000104284
190	0.000103724
191	0.000103228
192	0.00010267
193	0.000102172
194	0.000101614
195	0.000101118
196	0.00010062
197	0.000100062
198	0.000099566
199	0.000099006
200	0.00009851
201	0.000098014
202	0.000097454
203	0.00009702
204	0.000096462
205	0.000095966
206	0.000095468
207	0.00009491
208	0.000094414
209	0.000093916
210	0.00009342
211	0.000092924
212	0.000092428
213	0.000091868
214	0.000091372
215	0.000090876
216	0.000090378
217	0.000089882
218	0.000089386
219	0.000088888
220	0.000088392
221	0.000087958
222	0.0000874
223	0.000086902
224	0.000086406
225	0.00008591
226	0.000085474
227	0.00008504
228	0.000084544
229	0.000084048
230	0.00008355
231	0.000083054
232	0.000082558
233	0.00008206
234	0.000081626
235	0.00008113
236	0.000080634
237	0.000080198
238	0.000079702
239	0.000079268
240	0.00007877
241	0.000078336
242	0.00007784
243	0.000077344
244	0.000076908
245	0.000076474
246	0.00007604
247	0.000075544
248	0.000075046
249	0.000074612
250	0.000074178
251	0.00007368
252	0.000073246
253	0.000072812
254	0.000072316
255	0.00007188
256	0.000071446
257	0.000071012
258	0.000070516
259	0.00007008
260	0.000069646
261	0.000069212
262	0.000068778
263	0.000068342
264	0.000067908
265	0.000067474
266	0.00006704
267	0.000066604
268	0.00006617
269	0.000065736
270	0.000065302
271	0.000064866
272	0.000064494
273	0.000064122
274	0.000063626
275	0.000063252
276	0.000062818
277	0.000062384
278	0.00006195
279	0.000061514
280	0.00006108
281	0.000060708
282	0.000060274
283	0.000059838
284	0.000059466
285	0.000059032
286	0.00005866
287	0.000058224
288	0.000057852
289	0.000057418
290	0.000057046
291	0.000056672
292	0.000056238
293	0.000055866
294	0.000055494
295	0.000055122
296	0.000054686
297	0.000054314
298	0.000053942
299	0.00005357
300	0.000053196
301	0.000052762
302	0.00005239
303	0.000052018
304	0.000051644
305	0.000051272
306	0.0000509
307	0.000050528
308	0.000050156
309	0.000049782
310	0.000049472
311	0.0000491
312	0.000048728
313	0.000048356
314	0.000047982
315	0.00004761
316	0.0000473
317	0.000046928
318	0.000046554
319	0.000046182
320	0.000045872
321	0.0000455
322	0.00004519
323	0.000044816
324	0.000044506
325	0.000044134
326	0.000043824
327	0.000043514
328	0.00004314
329	0.00004283
330	0.000042458
331	0.000042148
332	0.000041838
333	0.000041528
334	0.000041154
335	0.000040844
336	0.000040534
337	0.000040224
338	0.000039914
339	0.000039602
340	0.000039292
341	0.000038982
342	0.000038672
343	0.000038362
344	0.00003805
345	0.00003774
346	0.00003743
347	0.00003712
348	0.00003681
349	0.000036562
350	0.00003625
351	0.00003594
352	0.000035692
353	0.000035382
354	0.000035072
355	0.000034824
356	0.000034512
357	0.000034264
358	0.000033954
359	0.000033706
360	0.000033458
361	0.000033148
362	0.000032898
363	0.000032588
364	0.00003234
365	0.000032092
366	0.000031782
367	0.000031534
368	0.000031284
369	0.000031036
370	0.000030788
371	0.000030478
372	0.00003023
373	0.000029982
374	0.000029734
375	0.000029484
376	0.000029236
377	0.000028988
378	0.00002874
379	0.000028492
380	0.000028306
381	0.000028058
382	0.000027808
383	0.00002756
384	0.000027312
385	0.000027126
386	0.000026878
387	0.00002663
388	0.000026444
389	0.000026194
390	0.000025946
391	0.00002576
392	0.000025512
393	0.000025326
394	0.00002514
395	0.000024892
396	0.000024706
397	0.000024456
398	0.00002427
399	0.000024084
400	0.000023898
401	0.00002365
402	0.000023464
403	0.000023278
404	0.000023092
405	0.000022906
406	0.000022656
407	0.00002247
408	0.000022284
409	0.000022098
410	0.000021912
411	0.000021726
412	0.00002154
413	0.000021354
414	0.000021166
415	0.000021042
416	0.000020856
417	0.00002067
418	0.000020484
419	0.000020298
420	0.000020174
421	0.000019988
422	0.000019802
423	0.000019678
424	0.00001949
425	0.000019304
426	0.00001918
427	0.000018994
428	0.00001887
429	0.000018684
430	0.00001856
431	0.000018374
432	0.00001825
433	0.000018126
434	0.00001794
435	0.000017816
436	0.000017628
437	0.000017504
438	0.00001738
439	0.000017256
440	0.00001707
441	0.000016946
442	0.000016822
443	0.000016698
444	0.000016574
445	0.00001645
446	0.000016264
447	0.00001614
448	0.000016014
449	0.00001589
450	0.000015766

OpenFC代码 ：

1. 
   
2. !using["XSLSF"];                //使用命名空间XSLSF
   
3. //函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。
   
4. //t为自变量，A,B,C,D,E,F为函数值，dA,dB,dC,dD,dE,dF为右端函数值（即微分方程的值）。
   
5. //该函数被调用时将用到参数k1, k2, k4, k11, k16, k17，这6个参数正好是拟合参数，用模块变量传递这6个参数。
   
6. f(t,A,B,C,D,E,F,dA,dB,dC,dD,dE,dF::k1, k2, k4, k11, k16, k17)=
   
7. {
   
8.   dA=56*k1*D-k2*A*B,
   
9.   dB=-k2*A*B,
   
10.   dC=k2*A*B-k4*C*D-k16*C*E-56*k17*C*F,
    
11.   dD=0.002*k11*A- 56*k1*D,
    
12.   dE=-k16*C*E,
    
13.   dF= k16*C*E-56*k17*C*F
    
14. };
    
15. //用连分式法对微分方程组进行积分，获得理论值。
    
16. //t1,t2为积分的起点和终点。
    
17. //h,s为自动变量。
    
18. //模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。
    
19. 积分(t1,t2:h,s:hf,Array,step,eps)=
    
20. {
    
21.     s=ceil[(t2-t1)/step],        //计算积分步数
    
22.     h=(t2-t1)/s,                 //重新计算积分步长
    
23.     {   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
    
24.         t1=t1+h                  //积分步长增加
    
25.     }.until[abs(t1-t2)<h/2]      //连续积分至t2
    
26. };
    
27. //目标函数定义，自变量_k,_b为需要优化的参数，需要将这些参数传递给对应的模块变量k,b。
    
28. //模块变量：Array为工作数组；数组tEArray存放实验数据(ti,Ei)；max为实验数据组数；k,b为优化变量。
    
29. 优化(_k1, _k2, _k4, _k11, _k16, _k17:i,s,E,a1,a2,a3,a4,a5,a6,t1,t2:Array,tEArray,max,k1, k2, k4, k11, k16, k17)=
    
30. {
    
31.   k1=_k1, k2=_k2, k4=_k4, k11=_k11, k16=_k16, k17=_k17,  //传递优化变量，函数f中要用到k1, k2, k4, k11, k16, k17
    
32.   Array.setra(0,0,1.25*10^(-3),0,10^(-4),2*10^(-4),0),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
33.   i=1,s=0,t1=0,
    
34.   (i<max).while{
    
35.     tEArray.getra(i*2,&t2,&E),
    
36.     积分(&t1,t2),                    //从t1积分到t2
    
37.     Array.getra(0,a1,a2,a3,a4,&a5,a6),  //从数组Array获得t2时的积分值
    
38.     s=s+[a5-E]^2,                    //累加理论值与实验值差的平方
    
39.     i++
    
40.   },
    
41.   s
    
42. };
    
43. 验证(_k1, _k2, _k4, _k11, _k16, _k17:i,s,E,a1,a2,a3,a4,a5,a6,t1,t2:Array,tEArray,max,k1, k2, k4, k11, k16, k17)=
    
44. {
    
45.   k1=_k1, k2=_k2, k4=_k4, k11=_k11, k16=_k16, k17=_k17,  //传递优化变量，函数f中要用到k1, k2, k4, k11, k16, k17
    
46.   printff{"\r\n               t             A模拟值             B模拟值          C模拟值          D模拟值          E模拟值          F模拟值          E实验值\r\n"},
    
47.   Array.setra(0,0,1.25*10^(-3),0,10^(-4),2*10^(-4),0),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
48.   i=1,s=0,t1=0,
    
49.   (i<max).while{
    
50.     tEArray.getra(i*2,&t2,&E),
    
51.     积分(&t1,t2),                    //从t1积分到t2
    
52.     Array.getra(0,&a1,&a2,&a3,&a4,&a5,&a6),            //从数组Array获得t2时的积分值
    
53.     s=s+[a5-E]^2,                    //累加理论值与实验值差的平方
    
54.     printff{"{1,r,18.8}{2,r,18.8}{3,r,18.8}{4,r,18.8}{5,r,18.8}{6,r,18.8}{7,r,18.8}{8,r,18.8}\r\n",t1,a1,a2,a3,a4,a5,a6,E},
    
55.     i++
    
56.   },
    
57.   printff{"\r\n目标函数终值={1,r}\r\n\r\n",s}
    
58. };
    
59. main(:x,xx,g,d,u,v,k,i,t1,t2,t3,t4,t5,t6,t7,_eps:Array,tEArray,max,hf,step,eps)=
    
60. {
    
61.   hf=HFor("f"),                                    //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
    
62.   step=0.05,eps=1e-6,                              //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
    
63.   Array=new[rtoi(real_s),rtoi(6*15)],              //申请工作数组
    
64.   max=451,                                         //实验数据组数
    
65.   tEArray=new{rtoi(real_s),rtoi(max),rtoi(2),rtoi(EndType), //存放实验数据(ti,Ei)
    
66. 0 , 0.0002 ,
    
67. 1 , 0.000199566 ,
    
68. 2 , 0.000199814 ,
    
69. 3 , 0.000199068 ,
    
70. ...... //数据略
    
71. 447 , 0.00001614 ,
    
72. 448 , 0.000016014 ,
    
73. 449 , 0.00001589 ,
    
74. 450 , 0.000015766
    
75.     },
    
76.   x=new[rtoi(real_s),rtoi(7)],                     //申请工作数组
    
77.   xx=new[rtoi(real_s),rtoi(6),rtoi(7)],            //申请工作数组
    
78.   g=new[rtoi(real_s),rtoi(7)],                     //申请工作数组
    
79.   _eps=1e-100, d=2,u=1.1,v=0.1,k=100,            //变换d、u、v进一步求解，k为允许的最大迭代次数
    
80.   i=jsim[HFor("优化"),d,u,v,x,_eps,k,xx,g],        //求n维极值的单形调优法
    
81.   x.getra(0,&t1,&t2,&t3,&t4,&t5,&t6,&t7),          //获得最优参数值及目标终值
    
82.   printff{"\r\n实际迭代次数={1,i}, k1={2,r}, k2={3,r}, k4={4,r}, k11={5,r}, k16={6,r}, k17={7,r}, 目标函数终值={8,r}\r\n",i,t1,t2,t3,t4,t5,t6,t7},
    
83.   验证(t1,t2,t3,t4,t5,t6),                         //验证计算结果
    
84.   delete[Array],delete[tEArray],delete[x],delete[xx],delete[g]
    
85. };
    

复制代码

部分结果：

1. 
   
2. 实际迭代次数=100, k1=11.572119058286493, k2=20.605987644786079, k4=5.6399027460131057, k11=10.145516619333584, k16=35.565808626819823, k17=9.9756597283905499, 目标函数终值=1.1440703884564055e-009
   
3.                t             A模拟值             B模拟值          C模拟值          D模拟值          E模拟值          F模拟值          E实验值
   
4.                 1.    9.9455987e-005    1.2474379e-003    2.5529959e-006    3.0729319e-009    1.9999091e-004    9.0845252e-009      1.99566e-004
   
5.                 2.    9.8921631e-005    1.2448909e-003    5.0728026e-006    3.0616921e-009    1.9996377e-004    3.6176865e-008      1.99814e-004
   
6.                 3.    9.8395292e-005    1.2423627e-003    7.5558911e-006    3.0499714e-009    1.9991885e-004    8.0893664e-008      1.99068e-004
   
7.                 4.    9.7876749e-005     1.239853e-003    1.0002646e-005     3.123494e-009    1.9985641e-004    1.4277954e-007      1.99006e-004
   
8.                 5.    9.7366051e-005    1.2373614e-003    1.2413381e-005     3.099665e-009    1.9977674e-004    2.2131168e-007      1.98448e-004
   
9.                 6.    9.6862995e-005    1.2348878e-003    1.4788347e-005    3.0551383e-009    1.9968011e-004    3.1590477e-007      1.98634e-004
   
10.                 7.    9.6367467e-005    1.2324318e-003    1.7127737e-005    2.9752671e-009    1.9956679e-004    4.2591622e-007      1.97704e-004
    
11. ... ...
    
12.                91.    7.3024114e-005    1.0689567e-003    1.0863441e-004    2.3030431e-009     1.585819e-004    1.0427352e-005      1.59032e-004
    
13.                92.    7.2898677e-005    1.0673508e-003    1.0899453e-004      2.29691e-009    1.5796936e-004    1.0406658e-005       1.5841e-004
    
14.                93.    7.2775858e-005      1.06575e-003    1.0934911e-004    2.2911529e-009    1.5735718e-004    1.0384815e-005      1.57852e-004
    
15.                94.    7.2655673e-005    1.0641543e-003    1.0969835e-004    2.2513477e-009    1.5674542e-004    1.0361888e-005      1.57232e-004
    
16.                95.    7.2538029e-005    1.0625636e-003    1.1004248e-004    2.2508112e-009    1.5613411e-004    1.0337937e-005      1.56672e-004
    
17.                96.     7.242289e-005    1.0609778e-003    1.1038171e-004    2.3065522e-009    1.5552329e-004    1.0313022e-005      1.56176e-004
    
18.                97.    7.2310356e-005    1.0593969e-003    1.1071626e-004    2.2978462e-009    1.5491301e-004    1.0287197e-005      1.55556e-004
    
19.                98.    7.2192075e-005    1.0578207e-003     1.110463e-004     1.056091e-008     1.543033e-004    1.0260515e-005      1.55122e-004
    
20.                99.    7.2092921e-005    1.0562493e-003    1.1137204e-004    2.2011907e-009     1.536942e-004    1.0233027e-005      1.54438e-004
    
21.               100.    7.1987866e-005    1.0546825e-003    1.1169367e-004     2.229601e-009    1.5308574e-004     1.020478e-005       1.5388e-004
    
22. ... ...
    

复制代码

TBE_Legend

第二个问题，将允许的最大迭代次数k＝500，运行了约40分钟，终于获得一组最优解：


实际迭代次数=213, k1=9.5111080627706936, k2=18.593155843279298, k4=-0.3629342128423555, k11=9.2785258647078628, k16= ...
wanglu 发表于 2009-11-28 13:23

maple得到的解析解

1. restart;
   
2. ode1:=diff(y(t),t)=k*y(t)*z(t)+0.095*b*z(t);
   
3. ode2:=diff(z(t),t)=-b*z(t)-0.222*z(t);
   
4. sol:=dsolve({ode1,ode2},{y(t),z(t)},build);
   
5. sol[1];sol[2];
   

复制代码

wanglu

maple得到的解析解
restart;
ode1:=diff(y(t),t)=k*y(t)*z(t)+0.095*b*z(t);
ode2:=diff(z(t),t)=-b*z(t)-0.222*z(t);
sol:=dsolve({ode1,ode2},{y(t),z(t)},build);
sol[1];sol[2];

TBE_Legend 发表于 2009-11-28 18:20

非常感谢TBE_Legend！

wanglu

第一个问题，太奇怪了，根据TBE_Legend给出的解析解，我带入初始值：
t＝0.01， y0＝2.329101563， z0＝5000
求得解析解中的参数c1和c2：
c1=-2.364794921874116
c2=5011.104980662317

然后修改代码：

1. !using["XSLSF"];                //使用命名空间XSLSF
   
2. //函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。
   
3. //t为自变量，y和z为函数值，dy和dz为右端函数值（即微分方程的值）。
   
4. //该函数被调用时将用到参数k和b，这2个参数正好是拟合参数，用模块变量传递这2个参数。
   
5. f(t,y,z,dy,dz::k,b)=
   
6. {
   
7.   dy=k*y*z+0.095*b*z,
   
8.   dz=-b*z-0.222*z
   
9. };
   
10. //用连分式法对微分方程组进行积分，获得理论值。
    
11. //t1,t2为积分的起点和终点。
    
12. //h,s为自动变量。
    
13. //模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。
    
14. 积分(t1,t2:h,s:hf,Array,step,eps)=
    
15. {
    
16.     s=ceil[(t2-t1)/step],        //计算积分步数
    
17.     h=(t2-t1)/s,                 //重新计算积分步长
    
18.     {   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步
    
19.         t1=t1+h                  //积分步长增加
    
20.     }.until[abs(t1-t2)<h/2]      //连续积分至t2
    
21. };
    
22. //目标函数定义，自变量_k,_b为需要优化的参数，需要将这些参数传递给对应的模块变量k,b。
    
23. //模块变量：Array为工作数组；数组tArray存放实验数据ti；数组yArray存放实验数据yi；max为数组tArray和yArray的长度；k,b为优化变量。
    
24. 优化(_k,_b:i,s,y,z,yy,t1,t2:Array,tArray,yArray,max,k,b)=
    
25. {
    
26.   k=_k,b=_b,                         //传递优化变量，函数f中要用到k,b
    
27.   Array.setra(0,2.329101563,5000),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
28.   i=1,s=0,t1=0.01,
    
29.   (i<max).while{
    
30.     tArray.getra(i,&t2), yArray.getra(i,&yy),
    
31.     积分(&t1,t2),                    //从t1积分到t2
    
32.     Array.getra(0,&y,z),             //从数组Array获得t2时的积分值
    
33.     s=s+[y-yy]^2,                    //累加理论值与实验值差的平方
    
34.     i++
    
35.   },
    
36.   s
    
37. };
    
38. y理论(t:c1,c2,tt,k,b)= c1=-2.364794921874116, c2=5011.104980662317, k=1.4083232846396376e-004, b=-1.4666335917093112e-004,
    
39.      tt=500*k*c2*exp[-1/500*(500*b+111)*t]/(500*b+111), {-19/200*[b*exp(tt)]/k+c1}*exp(-tt);
    
40. z理论(t:c2,b)= c2=5011.104980662317, b=-1.4666335917093112e-004, c2*exp[-1/500*(500*b+111)*t];
    
41. 验证(_k,_b:i,y,z,yy,t1,t2:Array,tArray,yArray,max,k,b)=
    
42. {
    
43.   k=_k,b=_b,
    
44.   printff{"\r\n               t            y实验值            y模拟值         y理论值         z模拟值         z理论值\r\n"},
    
45.   Array.setra(0,2.329101563,5000),   //设置积分初值，通过模块变量Array传递，Array是一个数组
    
46.   i=1,t1=0.01,
    
47.   (i<max).while{
    
48.     tArray.getra(i,&t2), yArray.getra(i,&yy),
    
49.     积分(&t1,t2),                    //从t1积分到t2
    
50.     Array.getra(0,&y,&z),            //从数组Array获得t2时的积分值
    
51.     printff{"{1,r,18.8}{2,r,18.8}{3,r,18.8}{4,r,18.8}{5,r,18.8}{6,r,18.8}\r\n",t1,yy,y,y理论(t1),z,z理论(t1)},
    
52.     i++
    
53.   }
    
54. };
    
55. main(:x,xx,g,d,u,v,k,i,t1,t2,t3,_eps:Array,tArray,yArray,max,hf,step,eps)=
    
56. {
    
57.   hf=HFor("f"),                                    //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄
    
58.   step=0.01,eps=1e-6,                              //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1
    
59.   Array=new[rtoi(real_s),rtoi(2*15)],              //申请工作数组
    
60.   max=18,                                          //实验数据组数
    
61.   tArray=new{rtoi(real_s),rtoi(max),rtoi(EndType), //存放实验数据ti
    
62.     0.01,
    
63.     0.68,
    
64.     1.1,
    
65.     1.63,
    
66.     2.07,
    
67.     2.67,
    
68.     3.09,
    
69.     3.64,
    
70.     4.65,
    
71.     5.1,
    
72.     5.58,
    
73.     6.11,
    
74.     6.63,
    
75.     7.06,
    
76.     7.62,
    
77.     8.66,
    
78.     9.04,
    
79.     9.63
    
80.     },
    
81.   yArray=new{rtoi(real_s),rtoi(max),rtoi(EndType), //存放实验数据yi
    
82.     2.329101563,
    
83.     3.851292783,
    
84.     4.50093936,
    
85.     6.749172247,
    
86.     9.112062872,
    
87.     9.691657366,
    
88.     11.16928013,
    
89.     10.91451823,
    
90.     16.44279204,
    
91.     18.29615885,
    
92.     21.63989025,
    
93.     25.78611421,
    
94.     26.34282633,
    
95.     26.50581287,
    
96.     27.63951823,
    
97.     35.02757626,
    
98.     35.5562035,
    
99.     36.10393415
    
100.     },
     
101.   x=new[rtoi(real_s),rtoi(3)],                     //申请工作数组
     
102.   xx=new[rtoi(real_s),rtoi(2),rtoi(3)],            //申请工作数组
     
103.   g=new[rtoi(real_s),rtoi(3)],                     //申请工作数组
     
104.   _eps=1e-100, d=550,u=1.6,v=0.4,k=500,            //变换d、u、v进一步求解，k为允许的最大迭代次数
     
105.   i=jsim[HFor("优化"),d,u,v,x,_eps,k,xx,g],        //求n维极值的单形调优法
     
106.   x.getra(0,&t1,&t2,&t3),                          //获得最优参数值及目标终值
     
107.   printff{"\r\n实际迭代次数={1,i}, k={2,r}, b={3,r}, 目标函数终值={4,r}\r\n",i,t1,t2,t3},
     
108.   验证(t1,t2),                                     //验证计算结果
     
109.   delete[Array],delete[tArray],delete[yArray],delete[x],delete[xx],delete[g]
     
110. };

复制代码

结果：

1. 
   
2. 实际迭代次数=118, k=1.4083028531001023e-004, b=-1.4651243025576023e-004, 目标函数终值=24.291246604109279
   
3.                t            y实验值            y模拟值         y理论值         z模拟值         z理论值
   
4.               0.68         3.8512928          3.556226   -5.4436048e-002         4309.3888         4309.3892
   
5.                1.1         4.5009394         4.5089213   -9.6698167e-002         3925.9863          3925.987
   
6.               1.63         6.7491722         5.9133543       -0.15899987          3490.466         3490.4669
   
7.               2.07         9.1120629         7.2439206       -0.21802501         3165.8441         3165.8451
   
8.               2.67         9.6916574          9.277603       -0.30824134         2771.2754         2771.2765
   
9.               3.09          11.16928         10.832718       -0.37722813         2524.7175         2524.7187
   
10.               3.64         10.914518         13.002374       -0.47347691         2234.7062         2234.7074
    
11.               4.65         16.442792         17.253471       -0.66206191         1786.1058         1786.1071
    
12.                5.1         18.296159         19.204669       -0.74862018         1616.4036         1616.4048
    
13.               5.58          21.63989         21.291253       -0.84118462         1453.1215         1453.1227
    
14.               6.11         25.786114          23.57465       -0.94248014         1291.9228          1291.924
    
15.               6.63         26.342826         25.768943         -1.039823         1151.1573         1151.1585
    
16.               7.06         26.505813         27.533738        -1.1181126         1046.4158         1046.4169
    
17.               7.62         27.639518         29.747617        -1.2163246         924.16271         924.16377
    
18.               8.66         35.027576         33.557029        -1.3853178         733.74424          733.7452
    
19.               9.04         35.556204          34.84099         -1.442277         674.42212         674.42303
    
20.               9.63         36.103934         36.714599        -1.5253943         591.67804          591.6789
    

复制代码

z模拟值   与  z理论值 很接近。但 y实验值  接近  y模拟值  ，与  y理论值  相差却如此大？难道所解c1不对？

大家帮忙解一下c1，谢谢！

wanglu

实际迭代次数=118, k=1.4083028531001023e-004, b=-1.4651243025576023e-004, 目标函数终值=24.291246604109279

               t            y实验值            y模拟值         y理论值         z模拟值         z理论值
              0.68         3.8512928          3.556226   -5.4436048e-002         4309.3888         4309.3892
               1.1         4.5009394         4.5089213   -9.6698167e-002         3925.9863          3925.987
              1.63         6.7491722         5.9133543       -0.15899987          3490.466         3490.4669
              2.07         9.1120629         7.2439206       -0.21802501         3165.8441         3165.8451
              2.67         9.6916574          9.277603       -0.30824134         2771.2754         2771.2765
              3.09          11.16928         10.832718       -0.37722813         2524.7175         2524.7187
              3.64         10.914518         13.002374       -0.47347691         2234.7062         2234.7074
              4.65         16.442792         17.253471       -0.66206191         1786.1058         1786.1071
               5.1         18.296159         19.204669       -0.74862018         1616.4036         1616.4048
              5.58          21.63989         21.291253       -0.84118462         1453.1215         1453.1227
              6.11         25.786114          23.57465       -0.94248014         1291.9228          1291.924
              6.63         26.342826         25.768943         -1.039823         1151.1573         1151.1585
              7.06         26.505813         27.533738        -1.1181126         1046.4158         1046.4169
              7.62         27.639518         29.747617        -1.2163246         924.16271         924.16377
              8.66         35.027576         33.557029        -1.3853178         733.74424          733.7452
              9.04         35.556204          34.84099         -1.442277         674.42212         674.42303
              9.63         36.103934         36.714599        -1.5253943         591.67804          591.6789

wanglu

补充：求解c1和c2时，k=1.4083232846396376e-004, b=-1.4666335917093112e-004，取优化解。

wanglu

似乎用这组解求c1和c2，效果更好：
k=1.4083028531001023e-004, b=-1.4651243025576023e-004

wanglu

是我自己解错了，求解c1时忘了y0＝2.329101563，重新求得c1＝53.3071289062511

仅给出结果，还是比较理想的。

实际迭代次数=118, k=1.4083028531001023e-004, b=-1.4651243025576023e-004, 目标函数终值=24.291246604109279

               t            y实验值            y模拟值         y理论值         z模拟值          z理论值
              0.68         3.8512928          3.556226         3.5562263      4309.3887631      4309.3887566
               1.1         4.5009394         4.5089213         4.5089218      3925.9863184      3925.9863125
              1.63         6.7491722         5.9133543         5.9133549      3490.4660381      3490.4660328
              2.07         9.1120629         7.2439206         7.2439213      3165.8441308      3165.8441261
              2.67         9.6916574          9.277603         9.2776039      2771.2753669      2771.2753627
              3.09          11.16928         10.832718         10.832719      2524.7174885      2524.7174847
              3.64         10.914518         13.002374         13.002375      2234.7061716      2234.7061682
              4.65         16.442792         17.253471         17.253473      1786.1058194      1786.1058168
               5.1         18.296159         19.204669          19.20467      1616.4035666      1616.4035641
              5.58          21.63989         21.291253         21.291255      1453.1214789      1453.1214767
              6.11         25.786114          23.57465         23.574653      1291.9227832      1291.9227813
              6.63         26.342826         25.768943         25.768945      1151.1573467      1151.1573449
              7.06         26.505813         27.533738          27.53374      1046.4158252      1046.4158236
              7.62         27.639518         29.747617          29.74762      924.16271166      924.16271027
              8.66         35.027576         33.557029         33.557032      733.74424214      733.74424103
              9.04         35.556204          34.84099         34.840993      674.42211527      674.42211425
              9.63         36.103934         36.714599         36.714603       591.6780432       591.6780423

TBE_Legend

我很喜欢做这样的题目，第一个很容易，即使不用解析解，也能很容易地求出来。我比较关心的是第二个，mmtc算了半天，我给关了，太慢了。

但是我对你的结果也表示怀疑，因为在优化的过程中，要解ode，很可能遇到ode在某个时间点以后的解就不存在了（优化过程中mmtc的结果有这样的提示）。我对mmtc的ode的算法还是比较有信心的。

所以，我希望能够找一个简单点的（或有准确解答的），但是题目类型完全一样。或许可以把这个方程的一些系数改改什么的。

wanglu

我很喜欢做这样的题目，第一个很容易，即使不用解析解，也能很容易地求出来。我比较关心的是第二个，mmtc算了半天，我给关了，太慢了。

但是我对你的结果也表示怀疑，因为在优化的过程中，要解ode，很可能遇到ode ...
TBE_Legend 发表于 2009-11-29 16:20

第二个运算量确实大了，Forcal运行约40分钟，才给出了一组解，我估计也是局部最优解。我用的“求n维极值的单形调优法”需要尝试求解，但运行时间长就难以尝试了。

如果对结果表示怀疑，可以验证一下：用我得出的最优解带入微分方程，然后求模拟值，并与实验值比较，给出实验值与模拟值差的平方和，我就是想请大家验证一下这个结果，若全部给出验证数据，就更好了。

k1=9.5111080627706936, k2=18.593155843279298, k4=-0.3629342128423555, k11=9.2785258647078628, k16=42.278370674583023, k17=23.125269056444612

我还在想，优化结果不理想，是不是也有可能是公式模型不太合理？例如下面这个例子：http://www.ilovematlab.cn/viewthread.php?tid=50048

数学模型：
dA/dt=-k1*A*B+k_1*C*D
dB/dt=-k1*A*B+k-1*C*D-k2*C*B+k-2*E*D-k3*E*B+k-3*F*D
dC/dt=k1*A*B-k_1*C*D-k2*C*B+k_2*E*D
dD/dt=k1*A*B-k_1*C*D+k2*C*B-k_2*E*D+k3*E*B-k_3*F*D
dE/dt=k2*C*B-k_2*E*D-k3*E*B+k_3*F*D
dF/dt=k3*E*B-k-3*F*D
拟合参数：k1,k2,k3,k_1,k_2,k_3
实验数据（ti; Ai, Bi, Ci, Di, Ei, Fi）
0 0.3999 15.5780 0 0 0 0
5.3833 0.0346 11.2058 0.1986 4.3713 0.0817 0.0850
9.1833 0.0249 11.1680 0.1597 4.4097 0.1607 0.0546
12.9833 0.0261 11.0306 0.1198 4.5472 0.1639 0.0901
16.7667 0.0366 10.6350 0.0322 4.9419 1.1042 0.2269
20.5667 0.0329 10.5443 0.0350 5.0329 0.0643 0.2677
23.4167 0.0431 10.6538 0.0533 4.9241 0.0565 0.2470

此例Forcal拟合结果很不好，若有兴趣可以试一下？

有准确解答的例子似乎不好找，都是自己解决不了的问题才放到网上来。

能改方程的系数也是个好办法，似乎最终就是找一个最精确的公式吗。

wanglu

上面的例子在本论坛也有讨论：http://forum.simwe.com/thread-884975-1-2.html

napler821127

微分方程组参数拟合问题在Matlab里如何求解？Forcal是否能求解更为复杂的微分方程组？比如：
f1(t)/dt=A1(f2(t)/dt,t)
f2(t)/dt=A2(f1(t)/t,f2(t),t)
...
即在微分方程组内含有变量的导数，对于这种方程组的参数拟合如何进行？

wanglu

微分方程组参数拟合问题在Matlab里如何求解？Forcal是否能求解更为复杂的微分方程组？比如：
f1(t)/dt=A1(f2(t)/dt,t)
f2(t)/dt=A2(f1(t)/t,f2(t),t)
...
即在微分方程组内含有变量的导数，对于这种方程组的参数 ...
napler821127 发表于 2009-12-2 10:26

对应这种情况，就再嵌套一层解非线性方程组即可。
把f2(t)/dt和f2(t)/dt当成求解变量。
如果能得到f2(t)/dt和f2(t)/dt的解析表达式就不用解方程组了。

napler821127

Forcal程序对微分方程组参数拟合的理论基础是什么？能否给几篇参考文献看看？谢了？
如果想用MATLAB进行参数拟合该如何处理呢？

wanglu

Forcal程序对微分方程组参数拟合的理论基础是什么？能否给几篇参考文献看看？谢了？
如果想用MATLAB进行参数拟合该如何处理呢？
napler821127 发表于 2009-12-2 15:36

我还真没有考虑过理论基础这个问题，你看我这样回答能否满足你的要求：
《C常用算法程序集》第二版，徐士良主编，清华大学出版社
这本书中有很多算法，我用Forcal封装了，算法怎样实现的，有些我不太懂，但接口我能看明白，并已转换成了Forcal函数。
“微分方程组参数拟合”我是这样实现的：拟合函数采用“求n维极值的单形调优法”，其中目标函数定义中要调用“积分一步的连分式法”（也可以采用其他积分方法）对微分方程进行积分，每次积分前，拟合参数通过Forcal的模块变量传递给积分函数。如果“微分方程组内含有变量的导数”，我就考虑在微分方程组定义中解解非线性方程组（把f2(t)/dt和f2(t)/dt当成求解变量），如果解方程组时还有复杂计算，依次类推嵌套实现即可。Forcal函数没有任何限制，不会阻碍你使用任何东西。程序虽越来越复杂，但不会显著降低Forcal的效率。
至于“求n维极值的单形调优法”等算法怎样实现的，因我没有《C常用算法程序集》电子版，故不好给你提供资料。但可以从以下站点下载到源代码：http://xoomer.virgilio.it/forcal/xiazai/forcal9/forcal9code.rar
你可以找一下《C常用算法程序集》这本书，有详细的说明。
我不知道MATLAB如何解此类问题，等待MATLAB高手解答。

messenger

Matlab求解这种问题至少有4种方法

1、用ODE求解器+优化工具箱。但这种方法并不完善，例如像TBE_Legend说前面说的那样，ODE方程无解。

2、用Simulink Design Optimization工具箱。

3、用TOMLAB/PROPT工具箱。

4、用System Identification工具箱。

后三种方法的算法与第1种方法的算法本质上应该相同，只是集成了一下，所以其各部份的协同工作效果较好。

我还真没有考虑过理论基础这个问题，你看我这样回答能否满足你的要求：
《C常用算法程序集》第二版，徐士良主编，清华大学出版社
这本书中有很多算法，我用Forcal封装了，算法怎样实现的，有些我不太懂，但接口我能看明白，并已转换成了Forcal函数。
“微分方程组参数拟合”我是这样实现的：拟合函数采用“求n维极值的单形调优法”，其中目标函数定义中要调用“积分一步的连分式法”（也可以采用其他积分方法）对微分方程进行积分，每次积分前，拟合参数通过Forcal的模块变量传递给积分函数。如果“微分方程组内含有变量的导数”，我就考虑在微分方程组定义中解解非线性方程组（把f2(t)/dt和f2(t)/dt当成求解变量），如果解方程组时还有复杂计算，依次类推嵌套实现即可。Forcal函数没有任何限制，不会阻碍你使用任何东西。程序虽越来越复杂，但不会显著降低Forcal的效率。
至于“求n维极值的单形调优法”等算法怎样实现的，因我没有《C常用算法程序集》电子版，故不好给你提供资料。但可以从以下站点下载到源代码：http://xoomer.virgilio.it/forcal/xiazai/forcal9/forcal9code.rar
你可以找一下《C常用算法程序集》这本书，有详细的说明。
我不知道MATLAB如何解此类问题，等待MATLAB高手解答。
wanglu 发表于 2009-12-2 16:32

TBE_Legend

Matlab求解这种问题至少有4种方法

1、用ODE求解器+优化工具箱。但这种方法并不完善，例如像TBE_Legend说前面说的那样，ODE方程无解。
dsfsf
2、用Simulink Design Optimization工具箱。

3、用TOMLAB/PROPT工具 ...
messenger 发表于 2009-12-2 20:14

我不会matlab，只用过一点点，但我感觉matlab解决这两位函数调用函数的问题不是很好，反正我是不会，我也知道思路是rockwoods写的那几个积分方程求解思路，希望有matlab老用户能给出这个问题的解决方法，我在用mmtc解，但目前效果不理想，我已经向FreddyMusic求助，希望他能帮忙找到问题的解决办法，到时我会贴出结果来。

simulink作函数拟合本版也有人发过一个，好像，忘记了。不知道哪位知道连接？

最后两个工具箱没听说过。

wanglu

感谢messenger 和TBE_Legend 的说明，期待matlab和mmtc的结果。