修正过的一维FDTD等离子体MATLAB代码(公式修正)

格式：doc
大小：28.00 KB
文档页数：3

clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%系统参数

TimeT=3000;%迭代次数

KE=2000;%网格数目

kc=400;%源的位置

kpstart=600;%等离子体开始位置

kpstop=1000;%等离子体终止位置

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%物理参数

f=1.356e7;%电磁场频率

fpe=1e8;%等离子体频率

c0=3e8;%真空中波速

zdelta=1e-1;%网格大小

dt=zdelta/(2*c0);%时间间隔

u0=2e10;%碰撞频率

wpe=2*pi*fpe;%等离子体圆频率

epsz=1/(4*pi*9*10^9); % 真空介电常数

mu=1/(c0^2*epsz);%磁常数

ex_low_m1=0;

ex_low_m2=0;

ex_high_m1=0;

ex_high_m2=0;

a0=2*u0/dt+(2/dt)^2;

a1=-8/(dt)^2;

a2=-2*u0/dt+(2/dt)^2;

b0=wpe^2+2*u0/dt+(2/dt)^2;

b1=2*wpe^2-8/(dt)^2;

b2=wpe^2-2*u0/dt+(2/dt)^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%初始化电磁场

Ex=zeros(1,KE);

Ex_Pre=zeros(1,KE);

Hy=zeros(1,KE);

Hy1=zeros(1,KE);

Hy_Pre=zeros(1,KE);

Dx=zeros(1,KE);

Dx_Pre=zeros(1,KE);

Dx=Ex;

Dx1=Ex;

Dx2=Ex;

Dx3=Ex;

Ex1=Ex;

Ex2=Ex;

Ex3=Ex;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%开始计算

for T=1:TimeT

%%%保存前一时间的电磁场

Ex_Pre=Ex;

Hy_Pre=Hy;

%%%%中间差分计算Dx

for i=2:KE

Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));

end

%%%%%%%%加入源

Dx(kc)=sin(2*pi*f*T*dt);

%%%计算等离子区域之前和之后的电场Ex

for i=1:kpstart-1

Ex(i)=Dx(i)/epsz;

end

for i=kpstop+1:KE

Ex(i)=Dx(i)/epsz;

end

Dx3=Dx2;

Dx2=Dx1;

Dx1=Dx;

% %等离子体区域的电场(修改过，添加了公式中的介电常数epsz)

for i=kpstart:kpstop

Ex(i)=(1/b0)*(a0*(Dx1(i)/epsz)+a1*(Dx2(i)/epsz)+a2*(Dx3(i)/epsz)-b1*Ex1(i)-b2*Ex2(i));

end

Ex2=Ex1;

Ex1=Ex;

Ex(1)=ex_low_m2;

ex_low_m2=ex_low_m1;

ex_low_m1=Ex(2);

Ex(KE)=ex_high_m2;

ex_high_m2=ex_high_m1;

ex_high_m1=Ex(KE-1);

%%%%%%%%%%%%%%%%%%计算磁场

for i=1:KE-1

Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));

end

subplot(2,1,1)

plot(Ex);

% pause(0.00001);

% plot([400,400],[-2e11,2e11],'-r');

% axis([0 8000 -20 20]);

% axis([0 2000 -2e11 2e11]);

grid on;

subplot(2,1,2)

plot(Hy);

% axis([0 8000 -1 1]);

% axis([0 2000 -6e8 6e8]);

grid on;

pause(0.0001);

end

matlab练习程序（修正指数曲线拟合）

对于⼀般的指数曲线如：y=a*e^(k*t)，可以先对两边求对数得到：log(y) = log(a)+k*t 这样的曲线，然后⽤最⼩⼆乘来计算系数。

但是对于修正指数曲线如：y=k+a*b^t 这样的函数，没法直接求对数然后⽤最⼩⼆乘，因为有⼀个常数项k，这⾥可以利⽤三和法来计算系数。

对于曲线 y = k+a*b^t，三和法计算各系数公式如下：其中：

其中：

m=n/3，这⾥的n就是所有数据的总个数。

这种⽅法不⽌可以计算形如y = k+a*b^t这样的修正指数曲线，还可以计算形如y=k*a^(b^t)这样的compertz曲线和形如y=1/(k+a*b^t)这样的logistic曲

线，因为compertz曲线可以通过求log得到修正指数曲线，⽽logistic曲线可以通过求倒数得到修正指数曲线。

matlab代码如下：

clear all;

close all;

clc;

a=45.73423;

b=0.4234;

k=34.34534;

t=1:30; %t为整数没问题

t=1:0.1:30; %t为⼩数有问题，只能求出k

y=k+a*b.^t;

plot(t,y);

m=floor(length(t)/3);

S1 = sum(y(1:m));

S2 = sum(y(m+1:2*m));

S3 = sum(y(2*m+1:3*m));

b = ((S3-S2)/(S2-S1))^(1.0/m);

a = (S2-S1)*(b-1)/(b*(b^m-1)^2);

k = (1.0/m)*(S1 - a*b*(b^m-1)/(b-1));

y1=k+a*b.^t;

hold on;

plot(t,y1,'ro');

newy = y - k;

%log(newy) = log(a) + log(b)*t;

X = [t' ones(length(t),1)];

Y = log(newy)';

C = inv(X'*X)*X'*Y;

多种频谱校正方法及matlab代码

多种频谱校正方法

采样间隔归一化成1T,采样长度为N.这样FFT离散谱线为0,1)iXiN（,相应的

频率分辨率2/(1/)NfN.设FFT离散谱线局部极高谱线为m(为了数学上简洁,

假定从0开始,注意在MATLAB环境下数组实际操作的是从1开始),记频偏量.我

们需要使用谱线m和与之相邻一条次高谱线，记这连续两条谱线中左边一条序号为M(当

次高谱线在m左侧时1Mm,反之Mm).

下面列出若干算法的计算公式

1.加矩形窗的精确谱校正[1]

iiiXUjV

1()sin()()cos()MMMMoptMMVVMUUMKUU



211cos()

sin()

cos()

sin()optMM

optMMKMZVUM

KMZVUM









21cos()cos()()ZMZMMmZZ

2.加矩形窗情形,采用解析单频模型的幅值比校正[1,2]

1||()||||M

MMXMmXX



3.加汉宁窗情形,采用解析单频模型的幅值比校正[1,2]

12||||()||||MM

MMXXMmXX



4.加矩形窗情形,采用解析单频模型的复比值校正[3]

1Re()M

MMXMmXX



5.加汉宁窗情形,采用解析单频模型的复比值校正[3]

12()MM

MMXXMmXX



6.加矩形窗情形,采用解析单频模型的复合复比值校正[3]1

1Re()MmMMXMmXX



1mRmmX

XX

,1

111mmLmmmmXX

XXXX



0.5)0.5)mLmR（（

7.加汉宁窗情形,采用解析单频模型的复合复比值校正[3]

12Re()MMmMMXXMmXX



12mmRmmXX

DFP算法及Matlab程序

传播优秀Word版文档，希望对您有帮助，可双击去除！

1 / 4 作业二用DFP算法求解1212221422)(minxxxxxxf，取Tx110，10010H。

一、求解：

TTTgHpgxxxxxg)2,4(,)2,4(,)42,422()(00002121

（1）求迭代点x1

令32040)()(2000pxf，得)(0的极小值点410，

所以得：.)4,3(,)5.0,1(,2,1,5.0,201001010001TTTTggyxxsgpxx

于是，由DFP修正公式有4138388410010000000000001syssyHyHyyHHHTTTT

下一个搜索方向为TgHp6,851111

（2）求迭代点x2

令5.5458)()(2111pxf，得)(1的极小值点541

于是得：TTgpxx)0,0(,2,421112，所以：8,)2,4(2fxxT，

因Hesse阵TGxG4222)(为正定阵，)(xf为严格凸函数，所以x为整体

极小点。

二、DFP算法迭代步骤如下：

（1）给定初始点0x，初始矩阵0H（通常取单位阵），计算0g，令k=0，给定控制误差。

（2）令kkkgHp。

（3）由精确一维搜索确定步长k，)(min)(0kkkkkpxfpxf

（4）令kkkkpxx1。

（5）若kg，则1kxx停；

否则令kkkxxs1 ，kkkggy1 。

（6）由DFP修正公式得1kH。令k=k+1，转步骤（2）传播优秀Word版文档，希望对您有帮助，可双击去除！

三、DFP算法matlab程序实现

最优化方法之修正牛顿法matlab源码(含黄金分割法寻找步长)

revisenewton.m

syms x1 x2 x3 xx;

% f = x1*x1 +x2*x2 -x1*x2 -10*x1 -4*x2 + 60 ;

% f = x1^2 + 2*x2^2 - 2*x1 *x2 -4*x1 ;

f = 100 * (x1^2 - x2^2) + (x1 -1 )^2 ;

hessen = jacobian(jacobian(f , [x1,x2]),[x1,x2]) ;

gradd = jacobian(f , [x1,x2]) ;

X0 = [0,0]' ;

B = gradd' ;

x1 = X0(1);

x2 = X0(2);

A = eval(gradd) ;

% while sqrt( A(1)^2 + A(2)^2) >0.1

i=0;

while norm(A) >0.1

i = i+1 ;

fprintf('the number of iterations is: %d\n', i)

if i>10

break;

end

B1 = inv(hessen)* B ;

B2= eval(B1);

% X1 = X0 - B2

% X0 = X1 ;

f1= x1 + xx * B2(1);

f2= x2 + xx* B2(2);

% ff = norm(BB) ?

syms x1 x2 ;

fT=[subs(gradd(1),x1,f1),subs(gradd(2),x2,f2)];

ff = sqrt((fT(1))^2+(fT(2))^2);

MinData = GoldData(ff,0,1,0.01);

x1 = X0(1);

x2 = X0(2);

x1 = x1 + MinData * B2(1)

x2 = x2 + MinData * B2(2)

A = eval(gradd)

End

GoldData.m

function MiniData = GoldData( f,x0,h0,eps)

下载提示

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

修正过的一维FDTD等离子体MATLAB代码(公式修正)

合集下载

matlab练习程序（修正指数曲线拟合）

多种频谱校正方法及matlab代码

DFP算法及Matlab程序

最优化方法之修正牛顿法matlab源码(含黄金分割法寻找步长)

文档推荐

最新文档