本帖最后由 eleanorexu 于 2013-3-28 20:49 编辑

我想求解的函数是：deltaL=n^2*h/sqrt(n^2-(sin(theta1))^2)-n^2*h/sqrt(n^2-(sin(theta2))^2

其中，deltaL，h，theta1，theta2都已知(都可测量)，想求n是多少，不知如何解。。。

此外最重要的是，还想求n对deltaL，h，theta1，theta2每一项的偏微分是多少，从而计算不确定度dn。

我的程序是：

********************************************************************************************

clear

clc

syms n deltaL theta1 theta2 h; %创建符号变量

n_exp = solve(deltaL-n^2*h/sqrt(n^2-(sin(theta1))^2)+n^2*h/sqrt(n^2-(sin(theta2))^2),n);  %计算出方程的根

%下面对每一个变量求偏导

d_h = diff(n_exp,h);

d_deltaL = diff(n_exp,deltaL);

d_theta1= diff(n_exp,theta1);

d_theta2= diff(n_exp,theta2);

%利用得到的表达式带入数值计算

%后缀为num的为真实数值

h_num = 0.01; %h(厚度)10mm

theta1_num = 57/180*pi; %光线和玻璃夹角初始值57度

theta2_num = 58/180*pi;%转了1度

deltaL_num = 9.322878675087701e-05;  %相应的光程变化9.322878675087701e-05m

%以下代表不确定度

dh = 1*10^-6;

ddeltaL=79*10^-9;

dtheta1 = 1/648000*pi;

dtheta2 = 1/648000*pi;

%在此检查一下n，看看靠谱的

n_num = subs(n_exp, {h,theta1,theta2,deltaL}, {h_num,theta1_num,theta2_num,deltaL_num});

%把偏导的数值算出来

d_h_num = subs(d_h, {h,theta1,theta2,deltaL}, {h_num,theta1_num,theta2_num,deltaL_num});

d_deltaL_num = subs(d_deltaL, {h,theta1,theta2,deltaL}, {h_num,theta1_num,theta2_num,deltaL_num});

d_theta1_num = subs(d_theta1, {h,theta1,theta2,deltaL}, {h_num,theta1_num,theta2_num,deltaL_num});

d_theta2_num = subs(d_theta2, {h,theta1,theta2,deltaL}, {h_num,theta1_num,theta2_num,deltaL_num});

%将偏导合起来得到不确定度数值

d_n_num =sqrt( (d_h_num .* dh).^2 + (d_deltaL_num .* ddeltaL).^2 + (d_theta1_num .* dtheta1).^2 + (d_theta2_num .* dtheta2).^2);

********************************************************************************************

但matlab没有给我返回一个具体的数字，命令窗口显示的是：

Warning: The solutions are parametrized by the symbols:

z1 = RootOf(z^5 + (z^4*(- h^2*sin(theta1)^2 + h^2*sin(theta2)^2 + deltaL^2))/(2*deltaL*h) -

(z^3*(2*deltaL*h*sin(theta1)^2 - 4*deltaL*h*sin(theta2)^2))/(2*deltaL*h) - (z^2*(deltaL^2*sin(theta1)^2 -

deltaL^2*sin(theta2)^2 - 2*h^2*sin(theta2)^4 + 2*h^2*sin(theta1)^2*sin(theta2)^2))/(2*deltaL*h) +

(z*(2*deltaL*h*sin(theta2)^4 - 2*deltaL*h*sin(theta1)^2*sin(theta2)^2))/(2*deltaL*h) + (h^2*sin(theta2)^6 -

h^2*sin(theta1)^2*sin(theta2)^4)/(2*deltaL*h), z) minus RootOf(2*deltaL*h*z^4 - h^2*z^3*sin(theta1)^2 +

h^2*z^3*sin(theta2)^2 + deltaL^2*z^3 + 3*deltaL*h*z^2*sin(theta2)^2 - 2*deltaL*h*z^2*sin(theta1)^2 -

h^2*z*sin(theta1)^2*sin(theta2)^2 - deltaL^2*z*sin(theta1)^2 + h^2*z*sin(theta2)^4 + deltaL^2*z*sin(theta2)^2 -

deltaL*h*sin(theta1)^2*sin(theta2)^2 + deltaL*h*sin(theta2)^4, z)

> In solve at 180

In Uncertainty at 5

******************************************************************************************

请问是因为所求的表达式太复杂了吗？theta1和theta2可以只差1秒，也就是说sin(theta1)和sin(theta2)可以相差很小，可以怎样化简再用matlab计算？

还是可以用原始公式算但有更好的计算方法？

谢谢了！！！

原始公式如下：  图片贴不上，给个附件吧。。。

2013-3-28 20:49 上传

点击文件名下载附件

3.49 KB, 下载次数: 1

包含原始公式