龙格库塔法
我写了一段代码是用龙格库塔法解微分方程的，程序总是报错，说是我在描述微分方程的表达式时出现了错误，请高手给看看%function main
fidin=fopen('lcl.txt');
fidout=fopen('dmatlab.txt','w');
while ~feof(fidin)
tline=fgetl(fidin);
if double(tline(1))>=48&&double(tline(1))<=57
fprintf(fidout,'%s\n\n',tline);
continue
end
end
fclose(fidout);
d=importdata('dmatlab.txt');
Tt=d(:,1);
U=d(:,2);
P=d(:,3);
Q=d(:,4);
NIND=40;%定义个体数
MAXGEN=100;%迭代代数
PRECI=34;%二进制位数
GGAP=0.9;%代沟
NVAR=4;%待辨识参数个数
N=249;
G=zeros(NIND,1);
B=zeros(NIND,1);
X=zeros(NIND,1);
TT=zeros(NIND,1);
Ps=zeros(NIND,1);
Qs=zeros(NIND,1);
E0=zeros(NIND,1);
C=zeros(NIND,1);
E=zeros(NIND,1);
E1=zeros(N,NIND);
x=zeros(N,1);
RR=zeros(NIND,1);
FieldD=,);rep(,);rep(,)];
Chrom=crtbp(NIND,NVAR*PRECI);%初始种群
gen=0;%代计数器
trace=zeros(MAXGEN,2);%迭代性能追踪
y=bs2rv(Chrom,FieldD);
G=y(:,1);
B=y(:,2);
X=y(:,3);
TT=y(:,4);
for i=1:NIND
Ps(i,1)=(U(1).^2)*G(i);
Qs(i,1)=(U(1).^2)*B(i);
end
Pd=P(1)-Ps;
Qd=Q(1)-Qs;
for i=1:NIND,
E0(i)=((Pd(i).^2+(Qd(i)-U(1)*U(1)/X(i)).^2).^0.5)*X(i)/U(1);
C(i)=E0(i)/((U(1)*U(1)-(Pd(i)*X(i)/E0(i)).^2).^0.5);
end
for j=1:NIND
f=inline('1/TT(j)*(-EE+C(j)*((U(1)*U(1)-(Pd(j)*X(j)/EE).^2).^0.5))','t','EE');%微分方程的右边项
%f=inline('t*EE','t','EE');
dt=0.004; %x方向步长
xleft=0.029000; %区域的左边界
xright=1.020991; %区域的右边界
xx=xleft:dt:xright; %一系列离散的点
n=length(xx); %点的个数
EE=E0(j);
%%(2)龙格库塔法
RK=EE;
for i=2:n
k1=f(xx(i-1),RK(i-1));
k2=f(xx(i-1)+dt/2,RK(i-1)+k1*dt/2);
k3=f(xx(i-1)+dt/2,RK(i-1)+k2*dt/2);
k4=f(xx(i-1)+dt,RK(i-1)+k3*dt);
RK(i)=RK(i-1)+dt*(k1+2*k2+2*k3+k4)/6;
RR(j)=RK(i);
end
end
运行后报的错误：：??? Error using ==> inlineeval
Error in inline expression ==> 1/TT(j)*(-EE+C(j)*((U(1)*U(1)-(Pd(j)*X(j)/EE).^2).^0.5))
??? Error using ==> eval
Undefined command/function 'TT'.
Error in ==> inline.subsref at 25
INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);
Error in ==> gyddj2 at 66
k1=f(xx(i-1),RK(i-1));
请高手指点
页:
[1]