三层介质瑞雷波程序

《三层介质瑞雷波程序》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、三层介质瑞雷波程序程序思想：1、首先根据凡友华的标量传递算法求出瑞雷波的频散函数，即令x5=0，它是关于波速和频率的非线性方程；2、然后固定频率的数值，使波速从1500到3500取值，给定某一搜索步长，搜索有根区间，只要找到一个有根区间，就调用二分法的求根子函数，求得零点；3、最后改变频率的数值重复2步，得到不同频率时候的零点值。主函数：[reslut.m]functionresult%求根的主函数i=1;%记录根的个数n=9.9;%n为搜索步长forf=0:10:3000;%频率分别在0到3000变化forct=1700:n:3500%波速在1700到3500之间变

2、化时搜素有根区间a=feval('rayleigh',ct,f);b=feval('rayleigh',ct+9.9,f);ifa==0nn(f/10+1,i)=ct;i=i+1;elseifa*b<0jie(f/10+1,i)=erfen('rayleigh',ct,ct+9.9,f);%调用求根函数i=i+1;p=i-1;elsecontinueendendi=1;endfid=fopen('test.txt','wt');forj=1:(3000/10)+1fprintf(fid,'%3.3e%3.3e',10*(j-1),jie(j,1:p));fprintf

3、(fid,'');endfclose(fid);二分法求根子函数：[erfen.m]functionc=erfen(fun,a,b,t)%二分法求解频散方程的子函数%传递参数时fun是一个函数文件%a,b区间两端点%delta计算精度%c所求方程的根delta=0.0001;ya=feval(fun,a,t);yb=feval(fun,b,t);max1=1+round((log(b-a)-log(delta))/log(2));%max1表示最大迭代次数fork=1:max1c=(b+a)/2;yc=feval(fun,c,t);ifyc==0a=c;b=c;e

4、lseifyb*yc>0b=c;yb=yc;elsea=c;ya=yc;endifb-a

5、1:n-1)=[5,2];%各层的厚度fork=1:nvp(k)=sqrt(ct^2/zb(k)^2-1);vs(k)=sqrt(ct^2/hb(k)^2-1);t(k)=1-ct^2/(2*hb(k)^2);g(k)=ct^2/(2*hb(k)^2);endfork=1:n-1r(k)=vp(k)^2;s(k)=vs(k)^2;p(k)=vp(k)*kk*h(k);q(k)=vs(k)*kk*h(k);a(k)=cos(p(k));b(k)=cos(q(k));c(k)=sin(p(k))/vp(k);d(k)=sin(q(k))/vs(k);l(k)=hb(k+1

6、)^2*md(k+1)/(hb(k)^2*md(k));endx1=1+vp(n)*vs(n);x2=t(n)+vp(n)*vs(n);x3=i*vp(n)*g(n);x4=-i*vp(n)*g(n);x5=-t(n)^2-vp(n)*vs(n);fork=2:-1:1;x1=x1/l(k);x5=x5*l(k);p1=-t(k)*x1+(1+t(k))*x2+x5;p2=x1-2*x2-x5;p3=g(k)*x3;p4=-g(k)*x4;p5=-t(k)^2*x1+2*t(k)*x2+x5;q1=a(k)*b(k)*p2+c(k)*d(k)*p5-a(k)*d(k)

7、*p3-b(k)*c(k)*p4;q2=a(k)*d(k)*s(k)*p2-b(k)*c(k)*p5+a(k)*b(k)*p3-c(k)*d(k)*s(k)*p4;q3=b(k)*c(k)*r(k)*p2-a(k)*d(k)*p5-c(k)*d(k)*r(k)*p3+a(k)*b(k)*p4;q4=c(k)*d(k)*r(k)*s(k)*p2+a(k)*b(k)*p5+b(k)*c(k)*r(k)*p3+a(k)*d(k)*s(k)*p4;x1=q1-q4+2*p1;x2=t(k)*q1-q4+(1+t(k))*p1;x3=g(k)*q2;x4=-g