關於部落格
  • 125505

    累積人氣

  • 20

    今日人氣

    0

    追蹤人氣

Octave/Matlab 一個數值分析程式比較完整的模式

 
核心-負責計算的程式, 是整個程式的心臟. 這裡指的是演算法本身, 寫好後是無須修改內容的.
   
問題-一個負責提供解題資訊的檔案, 對象是使用者. 內容包誇如:函數, 容忍值(TOL), 求解區間([a,b]), 初始值, ..., 等.
   
執行稿-呈現結果用的程式, 對象是客人. 這是一個經過包裝的介面, 呈現最後對方所需要的結果. 撰寫方式因人而異, 這是最能展現個人特色的部份.
[問題]
我們有很多的問題(習題). 我們把每個問題設計成個別的檔案, 供核心來讀取資訊. 以牛頓法為例, 一個問題需要有方程式 f(x)=0 的 f, f', TOL, initial gauss(初始值). 因此我們可以設計一個讓核心讀取的範例:
 
function <輸出元>=<問題名稱>(<控制元>,<變數, 如果有函數的話>)
%問題名稱 - 為問題命名.
%控制元 - 讓核心來控制他想要這個問題的哪個資訊. 這裡設定了 f, df, tol, p0 供核心選擇.
%變數 - 當問題要回傳函數值 f(x), f'(x) 時就需要 x 這個變數.

switch <控制元>
    case 'f'
       
<輸出元>=<f的方程式>;
    case 'df'
       
<輸出元>=<f'的方程式>;
    case 'tol'
       
<輸出元>=<容忍值>;
    case 'p0'
       
<輸出元>=<初始值>;
end
 
當問題設定好之後, 讀取的方式為
feval(@<問題名稱>, <控制元>, <變數>)
範例:(P65, Example1)
設定好問題後你就可以再命令列嘗試讀取各部份來確定問題可以正常輸出.
>> feval(@ex2301,'tol') %讀取TOL的值

ans =

  1.0000e-006

>> feval(@ex2301,'p0') %讀取p0的值

ans =

    0.7854

>> feval(@ex2301,'f',0) %讀取 f 在 0 點的值

ans =

     1

>> feval(@ex2301,'df',0) %讀取 f' 在 0 點的值

ans =

    -1
 
[核心]
設計好一個可供測試用的問題後, 就開始來製作核心. 以下以牛頓法為例:
function p=NT_MD(func)
% p=NT_MD(func) - 牛頓法
%   func - 問題資訊
%   p - 計算結果的輸出
% 緊接著 function 之下的註解區塊是說明檔案
% 可以在此告訴對方這個演算法的使用方式
% 使用者可以使用 help NT_MD 指令來瀏覽你的說明

N0=1e6;

% 使用 feval 從 func 中讀取計算所需的 tol, p0
TOL=feval(func,'tol');
p0=feval(func,'p0');

for iter=1:N0
    % 演算法主體, 使用
    % feval(func,'f',p0) 來存取 f(p0) 的值與
    % feval(func,'df',p0) 來存取 f'(p0) 的值與
    % 來計算所需的 p = p0 - f(p0)/f'(p0)
    p=p0-feval(func,'f',p0)/feval(func,'df',p0);
    
    % 當判別條件 |p-p0| < TOL 成立時, 解為 p
    % 並使用 return 指令結束這個 function
    if abs(p-p0)<TOL
        return
    end
    
    % 下一個步驟的準備, 在此將 p0 指定為本次所計算出來的 p 值
    p0=p;
end

% 這裡是當收斂失敗時輸出失敗訊息的部份
if iter==N0
    fprintf('Error after N0=%dn',N0);
end
 
如此一來你便可使用 p=NT_MD(@ex2301) 來用 p 存取問題 ex2301 的解. 執行看看:
>> p=NT_MD(@ex2301) %將 p 指定為計算解

p =

    0.7391

>> fprintf('The Solution p=%10.7fn',p); %將解得的 p 作格式化輸出
The Solution p= 0.7390851
 
[執行稿]
當所有的演算法都完成後, 就可以來設計執行稿(以下只是範例, 這個部份就可以發揮個人的創意, 甚至繪圖等都可以放在這個部份). 我的想法是將三個演算法的結果放在同一行好彼此作比較. 所以一題使用如下的區塊
% 以 p1~p3 存取三個演算法的結果
% 利用格式化輸出將結果放置成一排
p1=NT_MD(@ez2306a);
p2=secant(@ez2306a);
p3=falsepn(@ez2306a);
fprintf('%15.8f, %15.8f, %15.8fn',p1,p2,p3);
 
如此有以下結果
為了輸入方便, 在這個區塊前面加上一行 func=@ez2306a
% 指定 func 為問題 ez2306a
func=@ez2306a;

然後我希望每題前面可以顯示一下小題號, 修改一下fprintf最前面
% 使用 %3s 來放置小題號字串 '(a)'
fprintf('%3s, %15.8f, %15.8f, %15.8fn','(a)',p1,p2,p3);
修改完成的區塊整個如下:
 
% 指定 func 為問題 ez2306a
% 以 p1~p3 存取三個演算法的結果
% 利用格式化輸出將結果放置成一排
func=@ez2306a;
p1=NT_MD(func);
p2=secant(func);
p3=falsepn(func);
fprintf('%3s, %15.8f, %15.8f, %15.8fn','(a)',p1,p2,p3);
 
我們來檢視修改完成的結果:前面多了'(a)'
接著用一樣的方式將所有題目都擺上, 結果就像下面這樣
不過這樣誰知道哪個是哪個演算法的結果呢?所以上面在加上
% 標示每行所使用的演算法
fprintf('%3s, %15s, %15s, %15sn','Q','Newton','Secant','False Position');
 
所以最後的結果是:自己把檔案下載回去後切換到該資料夾執行 hw2s 吧XD
檔案下載:
其中NT_MD, secant, falsepn 為核心, ex~, ez~為問題, hw2s為執行稿.


相簿設定
標籤設定
相簿狀態