大家帮我求解下这个方程组。或者帮我看下程序也行。但是我觉得我用了全局变量,不太好。而且结果也不对
方程和参数都在附件里面
clear all
clc
t0=0;tN=50e-9;h=5e-12;t=t0:h:tN;Index=length(t);
Ax=zeros(1,Index);Ay=zeros(1,Index);Phx=zeros(1,Index);Phy=zeros(1,Index);N=zeros(1,Index);
n=zeros(1,Index);AxD=zeros(1,Index);AyD=zeros(1,Index);PhxD=zeros(1,Index);PhyD=zeros(1,Index);
Ax0=1e-9;Ay0=1e-6;N0=1e-6;n0=1e-6;Phx0=1e-6;Phy0=1e-6;AxD0=0;AyD0=0;PhxD0=0;PhyD0=0;
j=1;
global k a gammae gammas gammaa gammap tao w kx ky u
k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=0.5e9;gammap=1e9;tao=5e-9;w=2.2176e15;kx=10e9;ky=10e9;u=1.6;
m=1;
% for u=1:0.01:1.8
% u
for i=1:Index
t1=t0+h;
K11=function_Ax_sub(t0,Ax0,Ay0,N0,n0,Phx0,Phy0,AxD0,PhxD0);
K21=function_Ay_sub(t0,Ax0,Ay0,N0,n0,Phx0,Phy0,AyD0,PhyD0);
K31=function_N_sub(t0,Ax0,Ay0,N0,n0,Phx0,Phy0);
K41=function_n_sub_m(t0,Ax0,Ay0,N0,n0,Phx0,Phy0);
K51=function_Phx_sub(t0,Ax0,Ay0,N0,n0,Phx0,Phy0,AxD0,PhxD0);
K61=function_Phy_sub(t0,Ax0,Ay0,N0,n0,Phx0,Phy0,AyD0,PhyD0);
K12=function_Ax_sub(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2,AxD0,PhxD0);
K22=function_Ay_sub(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2,AyD0,PhyD0);
K32=function_N_sub(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2);
K42=function_n_sub_m(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2);
K52=function_Phx_sub(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2,AxD0,PhyD0);
K62=function_Phy_sub(t0+h/2,Ax0+h*K11/2,Ay0+h*K21/2,N0+h*K31/2,n0+h*K41/2,Phx0+h*K51/2,Phy0+h*K61/2,AyD0,PhyD0);
K13=function_Ax_sub(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2,AxD0,PhxD0);
K23=function_Ay_sub(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2,AyD0,PhyD0);
K33=function_N_sub(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2);
K43=function_n_sub_m(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2);
K53=function_Phx_sub(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2,AxD0,PhyD0);
K63=function_Phy_sub(t0+h/2,Ax0+h*K12/2,Ay0+h*K22/2,N0+h*K32/2,n0+h*K42/2,Phx0+h*K52/2,Phy0+h*K62/2,AyD0,PhyD0);
K14=function_Ax_sub(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63,AxD0,PhxD0);
K24=function_Ay_sub(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63,AyD0,PhyD0);
K34=function_N_sub(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63);
K44=function_n_sub_m(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63);
K54=function_Phx_sub(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63,AxD0,PhxD0);
K64=function_Phy_sub(t0+h,Ax0+h*K13,Ay0+h*K23,N0+h*K33,n0+h*K43,Phx0+h*K53,Phy0+h*K63,AyD0,PhyD0);
Ax1=Ax0+(h/6)*(K11+2*K12+2*K13+K14);
Ay1=Ay0+(h/6)*(K21+2*K22+2*K23+K24);
N1=N0+(h/6)*(K31+2*K32+2*K33+K34);
n1=n0+(h/6)*(K41+2*K42+2*K43+K44);
Phx1=Phx0+(h/6)*(K51+2*K52+2*K53+K54);
Phy1=Phy0+(h/6)*(K61+2*K62+2*K63+K64);
Ax(j)=Ax1;
Ay(j)=Ay1;
N(j)=N1;
n(j)=n1;
Phx(j)=Phx1;
Phy(j)=Phy1;
if(j*h<=tao)
AxD0=0;AyD0=0;PhxD0=0;PhyD0=0;
else
AxD0=Ax(j-floor(tao/h));AyD0=Ay(j-floor(tao/h));PhxD0=Phx(j-floor(tao/h));PhyD0=Phy(j-floor(tao/h));
end
Ax0=Ax1;Ay0=Ay1;N0=N1;n0=n1;Phx0=Phx1;Phy0=Phy1;
j=j+1;
t0=t1;
end
% Powerx(m)=mean(Ax.*conj(Ax));
% Powery(m)=mean(Ay.*conj(Ay));
% m=m+1
% end
% plot(u,Powerx,'-');
% hold on
% plot(u,Powery,':');
% hold on
%
x=mean(Ax.*conj(Ax))
y=mean(Ay.*conj(Ay))
% % %
subplot(2,1,1)
plot(t,Ax.*conj(Ax))
subplot(2,1,2)
plot(t,Ay.*conj(Ay))
function y =function_Ax_sub(t,Ax,Ay,N,n,Phx,Phy,AxD,PhxD)
%UNTITLED1 Summary of this function goes here
% Detailed explanation goes here
global k gammaa a kx w tao
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=(k*(N-1)-gammaa)*Ax-k*n*Ay*(sin(Phy-Phx)+a*cos(Phy-Phx))+kx*AxD*cos(Phx-PhxD+w*tao);
end
function y =function_Ay_sub(t,Ax,Ay,N,n,Phx,Phy,AyD,PhyD)
global k gammaa a ky w tao
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=(k*(N-1)+gammaa)*Ay-k*n*Ax*(sin(Phy-Phx)-a*cos(Phy-Phx))+ky*AyD*cos(Phy-PhyD+w*tao);
end
function y =function_N_sub(t,Ax,Ay,N,n,Phx,Phy)
%UNTITLED1 Summary of this function goes here
% Detailed explanation goes here
global gammae u
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=gammae*(-N*(1+Ax*conj(Ax)+Ay*conj(Ay))+u+2*n*Ax*Ay*sin(Phy-Phx));
end
function y =function_n_sub_m(t,Ax,Ay,N,n,Phx,Phy)
%UNTITLED1 Summary of this function goes here
% Detailed explanation goes here
global gammae gammas
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=-gammas*n-gammae*n*(Ax*conj(Ax)+Ay*conj(Ay))+2*gammae*N*Ax*Ay*sin(Phy-Phx);
end
function y =function_Phx_sub(t,Ax,Ay,N,n,Phx,Phy,AxD,PhxD)
global k gammap a kx w tao
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=k*a*(N-1)-gammap+k*n*Ay/Ax*(cos(Phy-Phx)-a*sin(Phy-Phx))-kx*AxD/Ax*sin(Phx-PhxD+w*tao);
end
function y =function_Phy_sub(t,Ax,Ay,N,n,Phx,Phy,AyD,PhyD)
%UNTITLED1 Summary of this function goes here
% Detailed explanation goes here
global k gammap a ky w tao
% k=300e9;a=3;gammae=1e9;gammas=50e9;gammaa=-0.1e9;gammap=6e9;tao=5e-9;w=2.2176e15;kx=0;ky=0;
y=k*a*(N-1)+gammap-k*n*Ax/Ay*(cos(Phy-Phx)+a*sin(Phy-Phx))-ky*AyD/Ay*sin(Phy-PhyD+w*tao);
end
2009-8-30 10:24 上传
点击文件名下载附件
34.5 KB, 下载次数: 495