adams_pmecme.m

微分方程数值解法源代码

function data = Adams_PMECME(x0,y0,h,N)
global diff_f;
%函数功能：四阶Adams预测－校正方法.
%输入参数：（x0，y0）－－初值；h－－步长；N－－计算次数.
%返 回 值：解曲线上的点，格式为二维列向量，第一列为自变量，第二列为函数值.
%调用关系：调用四阶龙格－库塔法函数Adams_RK_four_step(x0,y0,h,N)产生初值.

syms x y;
diff_f = input('请输入微分函数dy/dx:');

data(1,1) = x0;
data(1,2) = y0;

%四阶龙格－库塔法提供初值.
data_RK = Adams_RK_four_step(x0,y0,h,3);
data(2,:) = data_RK(2,:);
data(3,:) = data_RK(3,:);
data(4,:) = data_RK(4,:);
x_former = data(4,1);
y_former = data(4,2);
for i = 1:4
    x = data(i,1);
    y = data(i,2);
    dy(i) = eval(diff_f);
end

c0 = 0;
p0 = 0;
for i = 4:N
    x = x_former + h;
    p1 = y_former + h * (55 * dy(i) - 59 * dy(i-1) + 37 * dy(i-2) - 9 * dy(i-3)) / 24;
    y = p1 + 251 /270 * (c0 - p0);
    dm = eval(diff_f);
    c1 = y_former + h * (9 * dm + 19 * dy(i) - 5 * dy(i-1) + dy(i-2)) / 24;
    y = c1 - 19 / 270 * (c1 - p1);
    dy(i+1) = eval(diff_f);
    
    %保存数据.
    data(i+1,1) = x;
    data(i+1,2) = y;
    
    %更新数据
    x_former = x;
    y_former = y;
    c0 = c1;
    p0 = p1;
end