weifen_lgkt.m

可以方便的求解一阶微分方程

function yy=weifen_lgkt(A,B,x0,y0,fx)
% yy=weifen_lgkt(A,B,x0,y0,fx)
% 系统状态方程
% dy/dx-Ay=B*fx 微分方程,A,B为参数
% fx x的函数,可选参数,如不输入该参数, 则要在该函数体内编辑该函数
% x0,y0 状态变量初始值
%注：x0,y0可以为列向量,求解微分方程组
N=1000; %计算步数
h=0.002; %计算步长
x=x0;
y=y0;
%---------------------
for k=1:N
%求解微分方程
fx(:,k)=x0+(k-1)*h; %例子1，在此选f(x)=x 
k0=A*y+B*fx(:,k); %四阶龙格-库塔法则求解微分方程
k1=A*(y+h*k0/2)+B*fx(:,k);
k2=A*(y+h*k1/2)+B*fx(:,k);
k3=A*(y+h*k2)+B*fx(:,k);
y=y+(k0+2*k1+2*k2+k3)*h/6;
%--------------------------------------------------------------------------------------%
yy(k)=y; %输出
xx(k)=x0+k*h;

%y_jiexi(k)=2*fx(:,k)-2+4*exp(-fx(:,k)); %例子1的解析解
%y_e(k)=y_jiexi(k)-yy(k); %例子1解析解与数值解的误差
end
plot([x0 xx],[y0 yy])