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
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/278368.html
標籤：其他