MATLAB數值實驗：函數逼近法求方程的數值解

作者：凱魯嘎吉 - 博客園 http://www.cnblogs.com/kailugaji/

    這篇博客主要通過給定的數學迭代公式，利用MATLAB來迭代求解多項分數階微分方程的數值解，主要用到的是函數逼近法，一種是非線性化數值解法，一種為線性化數值解法，並繪制解析解與數值解的函數圖像，計算兩者的誤差。

1. 問題描述

2. MATLAB程序

demo_1.m

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

demo_2.m

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

3. 結果

>> demo_1
方法一中解析解與數值解之間的誤差為：0.169468
>> demo_2
方法二中解析解與數值解之間的誤差為：0.175177

方法一結果圖

方法二結果圖：

本博文問題來源：多項分數階常微分方程的數值微分法 http://max.book118.com/file_down/e37129207cb5d8d84fa03b346b308819.docx