MATLAB数值实验：函数逼近法求方程的数值解

clear
clc
format long % 数据形式为长精度
% Author: 凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/
%% 定义变量
alpha1 = 0.9;
alpha2 = 0.6;
alpha3 = 0.3; % 1>alpha1>alpha2>alpha3>0
%% 求解开始
T = 2; % 区间右端点
tau = 0.1; % 步长
TT = 0:tau:T; % t变量序列，也就是方程中的自变量t
N = length(TT)-1; % t变量序列个数-1
% 定义三个1*N的0矩阵用来储存方程中每一项的系数
b_alpha1 = zeros(1,N);
b_alpha2 = zeros(1,N);
b_alpha3 = zeros(1,N);
% 循环开始
for k = 0 : (N-1)
    b_alpha1(k+1) = ((1+k)^(1-alpha1)-(k)^(1-alpha1))/gamma(2-alpha1);
    b_alpha2(k+1) = ((1+k)^(1-alpha2)-(k)^(1-alpha2))/gamma(2-alpha2)*tau^(alpha1-alpha2);
    b_alpha3(k+1) = ((1+k)^(1-alpha3)-(k)^(1-alpha3))/gamma(2-alpha3)*tau^(alpha1-alpha3);
end
coe_0 = b_alpha1(0+1) + b_alpha2(0+1) + b_alpha3(0+1);
U = zeros(1,N+1); % 储存计算的结果
for n = 1:N
    temp = 0;
    for k = 0 : n-2
        temp = temp + (b_alpha1(n-k-1+1) + b_alpha2(n-k-1+1) + b_alpha3(n-k-1+1))*(U(k+1+1)-U(k+1));
    end
    temp0 = U(n);
    while 1
        temp1 = U(n-1+1) - temp /coe_0+ tau^(alpha1)*right_fun(TT(n+1),temp0,alpha1,alpha2,alpha3)/coe_0;
        % 计算误差 如果前一次迭代和后一次迭代的误差小于10^-7，那么久退出循环，并把最后一次迭代的值赋给U
        if abs(temp0-temp1) < 10^(-7)
            U(n+1) = temp1;
            break;
        else
            temp0 = temp1;
        end
    end
end
True_sol = true_fun(TT,alpha1); % 真实值
plot(TT,U,'-')
hold on
plot(TT,True_sol,'r*')
legend('数值解','解析解','Location','northwest')
title('Algorithm 1');
xlabel('t');
ylabel('u(t)');
err = max(abs(U-True_sol)); % 误差
saveas(gcf,sprintf('Algorithm 1.jpg'),'bmp'); %保存图片
fprintf('方法一中解析解与数值解之间的误差为：%f\n', err);
function aa = true_fun(t,alpha1)
aa = t.^(2+alpha1);
end
function bb = right_fun(t,u,alpha1,alpha2,alpha3)
bb = gamma(3+alpha1)/gamma(3)*t.^2+gamma(3+alpha1)/gamma(3+alpha1-alpha2)*t.^(2+alpha1-alpha2)+gamma(3+alpha1)/gamma(3+alpha1-alpha3)*t.^(2+alpha1-alpha3)+sin(t.^(2+alpha1))-sin(u);
end

clear
clc
format long
% Author: 凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/
%% 定义变量
alpha1 = 0.9;   
alpha2  = 0.6;  
alpha3 = 0.3;   % 1 > alpha1 > alpha2 > alpha3 > 0
%% 求解开始
T = 2;      
tau = 0.1;  
TT = 0:tau:T;
N = length(TT)-1;
% 定义三个1*N的0矩阵用来储存方程中每一项的系数
b_alpha1 = zeros(1,N); 
b_alpha2 = zeros(1,N); 
b_alpha3 = zeros(1,N);
% 循环开始
for k = 0 :(N-1)
    b_alpha1(k+1) = ((1+k)^(1-alpha1)-(k)^(1-alpha1))/gamma(2-alpha1);
    b_alpha2(k+1) = ((1+k)^(1-alpha2)-(k)^(1-alpha2))/gamma(2-alpha2)*tau^(alpha1-alpha2);
    b_alpha3(k+1) = ((1+k)^(1-alpha3)-(k)^(1-alpha3))/gamma(2-alpha3)*tau^(alpha1-alpha3);
end
coe_0 = b_alpha1(0+1) + b_alpha2(0+1) + b_alpha3(0+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%
U = zeros(1,N+1); 
temp0=U(2);
%% 第一个值特殊处理
while 1
    temp1= tau^(alpha1)*right_fun(TT(1+1),temp0,alpha1,alpha2,alpha3)/coe_0 ;
    if (abs(temp0-temp1) < 10^(-7) )
        U(2) = temp1;
        break;
    else
        temp0 = temp1;
    end
end
%% 2~N个值的计算
for n = 2:N
    temp = 0;
    for k = 0 : n-2
        temp = temp + (b_alpha1(n-k-1+1) + b_alpha2(n-k-1+1) + b_alpha3(n-k-1+1))*(U(k+1+1)-U(k+1));
    end
    temp0 = U(n);
    while 1
        XX=2*U(n-1+1)-U(n-2+1);
        temp1 = U(n-1+1) - temp /coe_0+ tau^(alpha1)*right_fun(TT(n+1),XX,alpha1,alpha2,alpha3)/coe_0;
        % 计算误差 如果前一次迭代和后一次迭代的误差小于10^-7，那么久退出循环，并把最后一次迭代的值赋给U
        if (abs(temp0-temp1) < 10^(-7) )
            U(n+1) = temp1;
            break;
        else
            temp0 = temp1;
        end
    end
end
True_sol = true_fun(TT,alpha1); % 真实值
plot(TT,U,'-')
hold on
plot(TT,True_sol,'r*')
legend('数值解','解析解','Location','northwest')
title('Algorithm 2');
xlabel('t');
ylabel('u(t)');
err = max(abs(U-True_sol)); % 误差
saveas(gcf,sprintf('Algorithm 2.jpg'),'bmp'); %保存图片
fprintf('方法二中解析解与数值解之间的误差为：%f\n', err);
function aa = true_fun(t,alpha1)
aa = t.^(2+alpha1);
end
function bb = right_fun(t,u,alpha1,alpha2,alpha3)
bb = gamma(3+alpha1)/gamma(3)*t.^2+gamma(3+alpha1)/gamma(3+alpha1-alpha2)*t.^(2+alpha1-alpha2)+gamma(3+alpha1)/gamma(3+alpha1-alpha3)*t.^(2+alpha1-alpha3)+sin(t.^(2+alpha1))-sin(u);
end

>> demo_1
方法一中解析解与数值解之间的误差为：0.169468
>> demo_2
方法二中解析解与数值解之间的误差为：0.175177