马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
一共八个方程,构成方程组。未知数是 x,c2,c3,c4,c5,c6,c7,c8。因为c1可以赋值为1.所以不用解。
我也是初学matlab,程序也是参考matlab的帮助写的,但是一直有问题,请大家帮忙看看。不胜感激。
% solve x c2 c3 c4 c5 c6 c7 c8
clear;
clc;
% =============== Define some parameters ================== %
Lambda = 1.55*1E-6;
Epsilon = 8.85*1E-12;
Mu = 4*pi*1E-7;
k = (2*pi)/Lambda;
r1 = 2*1E-5;
r2 = 3*1E-5;
n1 = 1.6;
n2 = 1.49;
syms x c2 c3 c4 c5 c6 c7 c8;
m=0;
Epsilon1 = Epsilon*n1^2;
Epsilon2 = Epsilon*n2^2;
u = (k^2*n1^2 - x^2)^0.5;
w = (x^2 - k^2*n2^2)^0.5;
f = 6*pi/Lambda*1E+8;
c1=1;
% ================= End of Definition ====================== %
% ================definition transform ===================== %
Dbslj = inline('1/2*(besselj(m-1,x)-besselj(m+1,x))');
Dbsly = inline('1/2*(bessely(m-1,x)-bessely(m+1,x))');
Dbsli = inline('1/2*(besseli(m-1,x)+besseli(m+1,x))');
Dbslk = inline('-1/2*(besselk(m-1,x)-besselk(m+1,x))');
% =============end of definition transform====================%
%================definition function==========================%
F1 = c1*besseli(m,w*r1)-c3*besselj(m,u*r1)-c4*bessely(m,u*r1);
F2 = c2*besseli(m,w*r1)-c5*besselj(m,u*r1)-c6*bessely(m,u*r1);
F3 = c3*besselj(m,u*r2)+c4*bessely(m,u*r2)-c7*besselk(m,w*r2);
F4 = c5*besselj(m,u*r2)+c6*bessely(m,u*r2)-c8*besselk(m,w*r2);
F5 = c1*(m*x)/(w^2*r1)*besseli(m,w*r1)+c2*(i/w)*Mu*f*Dbsli(m,w*r1)...
-c3*(m*x)/(u^2*r1)*besselj(m,u*r1)-c4*(m*x)/(u^2*r1)*bessely(m,u*r1)...
-c5*(i/u)*Mu*f*Dbslj(m,u*r1)-c6*(i/u)*Mu*f*Dbsly(m,u*r1);
F6 = -c1*(i*f*Epsilon2/w)*Dbsli(m,w*r1)+c2* (m*x)/(w^2*r1)*besseli(m,w*r1)...
+c3*(i*f*Epsilon1/u)*Dbslj(m,u*r1)+c4*(i*f*Epsilon1/u)*Dbsly(m,u*r1)...
-c5*(m*x)/(u^2*r1)*besselj(m,u*r1)-c6*(m*x)/(u^2*r1)*bessely(m,u*r1);
F7 = c3*(m*x)/(u^2*r2)*besselj(m,u*r2)+c4*(m*x)/(u^2*r2)*bessely(m,u*r2)...
+c5*(i*f*Mu/u)*Dbslj(m,u*r2)+c6*(i*f*Mu/u)*Dbsly(m,u*r2)...
-c7*(m*x)/(w^2*r2)*besselk(m,w*r2)-c8*(i*f*Mu/w)*Dbslk(m,w*r2);
F8 = -c3*(i*f*Epsilon1/u)*Dbslj(m,u*r2)-c4*(i*f*Epsilon1/u)*Dbsly(m,u*r2)...
+c5*(m*x)/(u^2*r2)*besselj(m,u*r2)+c6*(m*x)/(u^2*r2)*bessely(m,u*r2)...
+c7*(i*f*Epsilon2/w)*Dbslk(m,w*r2)-c8*(m*x)/(w^2*r2)*besselk(m,w*r2);
%==============================end of function======================%
%============================== program==========================%
function y = myFF(xc);
y=[F1;F2;F3;F4;F5;F6;F7;F8];
xc0 = input('[]');
xc = fsolve(@myFF,xc0,options)
如果把program写到一个m文件里,一直显示function定义有问题,我也尝试吧 program 新建一个m文件,但是解出来的结果根本不对。正确的结果是:x的值应该是介于6000000~6500000之间,有很多个值,也就是说,有很多组解。
求解方程.doc
(31 KB, 下载次数: 2)
|