[未答复] Adams预估校正法解分数阶常微分方程(等式右端函数线性可以,非线性为什么不可以?)

[复制链接]
catusH 发表于 2021-11-21 18:02:10
本帖最后由 catusH 于 2021-11-21 18:20 编辑

function Testpredictorcorrector
clc;
B=1;
alpha=0.5;
dt=0.01 ;   
t_max=1;  
t=0:dt:t_max;
n=length(t);
u=zeros(n,1);
u1=zeros(n,1);


y0=0;
b=zeros(n-1,1);
a=zeros(2*n-3,1);

for i=1:(n+1)
     T(i)=(i-1)/n*t_max;
%f(i)=gamma(6)*T(i)^(5-alpha)/gamma(6-alpha)-3*gamma(5)*T(i)^(3-alpha)/gamma(5-alpha)+gamma(5)*T(i)^(3-alpha)/gamma(4-alpha)+(T(i)^5-3*T(i)^4+2*T(i)^3)^2;
f(i)=T(i)^(3+alpha)+gamma(4+alpha)*T(i)^3/6;
end
f0=0;
for i=1:n
    b(i)=0;
     
    for j=0:(i-1)
        if j==0
          % b(i)=b(i)+(dt^alpha/alpha*((i-j)^alpha-(i-1-j)^alpha))*(f0-B*y0^2);
            b(i)=b(i)+(dt^alpha/alpha*((i-j)^alpha-(i-1-j)^alpha))*(f0-B*y0);% Prediction term
        else
           b(i)=b(i)+(dt^alpha/alpha*((i-j)^alpha-(i-1-j)^alpha))*(f(j)-B*y(j)^2);% Prediction term
           % b(i)=b(i)+(dt^alpha/alpha*((i-j)^alpha-(i-1-j)^alpha))*(f(j)-B*y(j));
        end
    end
    y(i)=y0+1/gamma(alpha)*b(i); % Prediction term
    a(i)=0;
    for j=0:i
        if j==0
           % a(i)=a(i)+dt^alpha/(alpha*(alpha+1))*((i-1)^(alpha+1)-(i-1-alpha)*(i)^alpha)*(f0-B*y0^2);% Correction term
            a(i)=a(i)+dt^alpha/(alpha*(alpha+1))*((i-1)^(alpha+1)-(i-1-alpha)*(i)^alpha)*(f0-B*y0);
        elseif j==i
            %a(i)=a(i)+dt^alpha/(alpha*(alpha+1))*(f(i)-B*y(i)^2); % Correction term
             a(i)=a(i)+dt^alpha/(alpha*(alpha+1))*(f(i)-B*y(i));
        else
            %a(i)=a(i)+(dt^alpha/(alpha*(alpha+1))*((i-1-j+2)^(alpha+1)+(i-1-j)^(alpha+1)-2*(i-1-j+1)^(alpha+1)))*(f(j)-B*y(j)^2); % Correction term
        a(i)=a(i)+(dt^alpha/(alpha*(alpha+1))*((i-1-j+2)^(alpha+1)+(i-1-j)^(alpha+1)-2*(i-1-j+1)^(alpha+1)))*(f(j)-B*y(j));
        end
      
    end
    y(i)=y0+1/gamma(alpha)*a(i); % Correction result
end
for i=n:(-1):1
    y(i+1)=y(i);
    %u1(i+1)=T(i+1)^5-3*T(i+1)^4+2*T(i+1)^3;%解析解
      u1(i+1)=T(i+1)^(3+alpha);
end
y(1)=y0;
u1(1)=y0;

   
%for i=1:n+1
  % if y(i)<0
   %  y(i)=0;
    %end
  % end


% hold on
  plot(T,y,'b-',T,u1,'--r')
hold on
%plot(T,y,'b-') % Plot the numerical result
title('线性函数 h=0.01')
% xlabel('t')
% ylabel('u(t)')
以上代码分别包含线性和非线性,主要在于a(i)项和b(i)项最后乘以的y的次数是一次还是高次
问题:
一次时也就是线性时正常计算,数值解与精确解吻合
非线性y^2时是数值不对的计算趋势
我的思路是在a(i)b(i)后面直接把y的次数平方了,这样是不是不对?




您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表