matlab 龙格库塔求解隐式方程,Matlab龙格库塔求解方程组问题

大家帮我求解下这个方程组。或者帮我看下程序也行。但是我觉得我用了全局变量，不太好。而且结果也不对

方程和参数都在附件里面

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