adams求解微分方程组

clear all;
clc;
x0 = 1;
b = 1; 
y0 = 1;
z0 = 1; 
t0 = 0;
h = 0.1;
df = @diff;
[k,X,Y] = Adams4t(df,t0,b,x0,y0,z0,h);
% temp = exact(X);
% plot(X,Y,'r',X,temp,'*');
plot(X,Y,'*');
legend('x','y','z');
% function result = exact(x)
%      result = 3./(1+x.^3);%exp(-5*x);
% end
% function result = diff(x,y)
%     result = -x.^2.*y.^2;%-5*x;%x.^3 - y./x;
% end
function result = diff(t,Y)
    x = Y(1);
    y = Y(2);
    z = Y(3);
    result(1) = y + 3*z +sin(5*t); 
    result(2) = x +cos(t);
    result(3) = x + z - 3*cos(3*t)*sin(4*t);
%     result(1) = -10*x + 10*y; 
%     result(2) = 28*x - y - x.*z;
%     result(3) = x.*y - 8/3*z;
end
function [k,T,Y] = Adams4t(funcfcn,t0,b,x0,y0,z0,h)
t = t0;
y = [x0,y0,z0];
p = 128; 
n = fix((b-t0)/h);
if n<5
	return;
end
T = zeros(p,1);%自变量是t，
Y = zeros(p,length(y)); f = zeros(p,length(y));
k = 1; T(k) = t0; Y(k,:) = y';
for k =2:4
	c1 = 1/6; c2 = 2/6;  c3 = 2/6;
	c4 = 1/6; a2 = 1/2;  a3 = 1/2;
	a4 = 1;   b21 = 1/2; b31 = 0;
	b32 = 1/2;b41 = 0;   b42 = 0;  b43 = 1;
	t1 = t+a2*h; t2 = t+a3*h;
	t3 = t+a4*h; 
    k1 = feval(funcfcn,t,y);
	y1 = y+b21*h*k1; 
    t = t+h;
	k2 = feval(funcfcn,t1,y1);
	y2 = y+b31*h*k1+b32*h*k2;
	k3 = feval(funcfcn,t2,y2);
    y3 = y+b41*h*k1+b42*h*k2+b43*h*k3;%y3 = y+h*k3;    
	k4 = feval(funcfcn,t3,y3);
	y = y+h*(c1*k1+c2*k2+c3*k3+c4*k4); %y = y+h*(1/6*k1+...+1/6*k4);
    T(k) = t; Y(k,:) = y';
end
for i = 1:4
    f(i,:) = feval(funcfcn,T(i),Y(i,:));
end
count_max = 100;%隐函数的最大迭代次数
 
for k = 4:n
    T(k+1) = T(1) + h*k;
    f(k+1,:) = feval(funcfcn,T(k+1),Y(k,:));
    count = 1;
    du = 1;
    
        while du>1e-6 && count <count_max
            for index = 1:length(y)
                s1(index) = Y(k,index) +(h/24)*sum(((f(k-2:k+1,index))'.*[1 -5 19 9]));
            end
            f(k+1,:) = feval(funcfcn,T(k+1),s1(:));
            for index = 1:length(y)
                s2(index) = Y(k,index) +(h/24)*sum(((f(k-2:k+1,index))'.*[1 -5 19 9]));
            end
            Y(k+1,:) = s2;
            du = max(abs(s2-s1));
            count = count+1;
        end
end
  %Y(k+1) = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
% for k = 4:n
% 	f(k) = feval(funcfcn,X(k),Y(k));
% 	X(k+1) = X(1) + h*k;
% 	Y(k+1) = Y(k) +(h/24)*sum(((f(k-3:k))'.*[-9 37 -59 55]));
%     
% % 	f(k+1) = feval(funcfcn,X(k+1),Y(k+1));%预估校正公式
% %     Y(k+1) = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
%  
% end
 
T = T(1:n+1);
Y = Y(1:n+1,:);
n = 1:n+1;
% wucha = wucha(1:n,:);
end

上一篇是求方程，这一篇是求方程组，采用的算例为四阶龙格库塔方法求解一次常微分方程组_HelloWorldTM的博客-CSDN博客，他是用RK，这里使用多步法