本帖最后由 daiqiangqin 于 2011-12-27 22:17 编辑

如题，用的这篇文章的公式编的M程序，程序和运行结构如下,x(i)为角度，y(i)为文章中的Z2.算出来的结果好像不符合，请大家给点意见。
不像某些人就想拿点东西出来在这里显摆，有种的就把程序和结果发出来给大家看看。
程序如下：
function varargout=saxplaxliu(varargin)
clc,clear
x0=0;xn=pi/4;y0=0;h=pi/160;
=lgkt4j(x0,xn,y0,h);
n=length(x);
fprintf('i x(i) y(i)\n');
for i=1:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i));
end
function z=f(x,y)
z=-((3.5+y)/(cos(x))^2)/((1.53419*cos(x)-sqrt(cos(x)))/(1.53419*sin(x)-sqrt(1-sin(x)))+tan(x));
function =lgkt4j(x0,xn,y0,h)
x=x0:h:xn;
n=length(x);
y1=x;
y1(1)=y0;
for i=1:n-1
K1=f(x(i),y1(i));
K2=f(x(i)+h/2,y1(i)+h/2*K1);
K3=f(x(i)+h/2,y1(i)+h/2*K2);
K4=f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
end
y=y1;
运行结果：
i x(i) y(i)
1 0.0000 0.0000
2 0.0196 0.1308
3 0.0393 0.2658
4 0.0589 0.4047
5 0.0785 0.5474
6 0.0982 0.6934
7 0.1178 0.8422
8 0.1374 0.9934
9 0.1571 1.1461
10 0.1767 1.2996
11 0.1963 1.4529
12 0.2160 1.6049
13 0.2356 1.7545
14 0.2553 1.9004
15 0.2749 2.0413
16 0.2945 2.1758
17 0.3142 2.3025
18 0.3338 2.4200
19 0.3534 2.5270
20 0.3731 2.6223
21 0.3927 2.7048
22 0.4123 2.7736
23 0.4320 2.8280
24 0.4516 2.8675
25 0.4712 2.8918
26 0.4909 2.9009
27 0.5105 2.8949
28 0.5301 2.8742
29 0.5498 2.8392
30 0.5694 2.7907
31 0.5890 2.7295
32 0.6087 2.6564
33 0.6283 2.5723
34 0.6480 2.4783
35 0.6676 2.3753
36 0.6872 2.2644
37 0.7069 2.1464
38 0.7265 2.0224
39 0.7461 1.8931
40 0.7658 1.7595
41 0.7854 1.6222
按夏的文章算出来的y(i)（也就是Z2）应该为负数。
大家找出原因可别忘了发上来，我算出来也会发上来。
ODE45算的结果：
x =
0
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0002
0.0004
0.0006
0.0008
0.0010
0.0019
0.0029
0.0038
0.0048
0.0096
0.0144
0.0192
0.0240
0.0436
0.0632
0.0829
0.1025
0.1221
0.1418
0.1614
0.1810
0.2007
0.2203
0.2399
0.2596
0.2792
0.2988
0.3185
0.3381
0.3577
0.3774
0.3970
0.4167
0.4363
0.4559
0.4756
0.4952
0.5148
0.5345
0.5541
0.5737
0.5934
0.6130
0.6326
0.6523
0.6719
0.6915
0.7112
0.7308
0.7445
0.7581
0.7718
0.7854
y =
0
0.0001
0.0001
0.0002
0.0002
0.0005
0.0007
0.0010
0.0012
0.0025
0.0037
0.0050
0.0062
0.0125
0.0188
0.0252
0.0315
0.0633
0.0953
0.1276
0.1601
0.2960
0.4358
0.5792
0.7259
0.8753
1.0269
1.1798
1.3333
1.4865
1.6381
1.7870
1.9319
2.0715
2.2044
2.3292
2.4444
2.5490
2.6415
2.7211
2.7868
2.8380
2.8742
2.8951
2.9009
2.8916
2.8677
2.8297
2.7783
2.7144
2.6388
2.5525
2.4564
2.3516
2.2390
2.1196
1.9944
1.9043
1.8121
1.7180
1.6222

• 近朗伯光型LED透镜的光学设计.pdf