romberg演算法matlab

『壹』 用龍貝格法求積分

先用另外2種方法。

format long
%【1】精確值。符號積分
it=int('(2/sqrt(pi))*exp(-x)',0,1)
Accurate=eval(it)
y=inline('(2/sqrt(pi))*exp(-x)')
%【2】Simpson方法
Simpson=quad(y,0,1)
delta=Simpson-Accurate

結果：
Accurate = 0.713271669674918

y = Inline function:
y(x) = (2/sqrt(pi))*exp(-x)

Simpson = 0.713271671228492

delta = 1.553574158208448e-009

【3】從網上找到一個，存為romberg.m
%=================
function R = romberg(f, a, b, n)
format long
% ROMBERG -- Compute Romberg table integral approximation.
%
% SYNOPSIS:
% R = romberg(f, a, b, n)
%
% DESCRIPTION:
% Computes the complete Romberg table approximation to the integral
%
% / b
% I = | f(x) dx
% / a
%
% PARAMETERS:
% f - The integrand. Assumed to be a function callable as
% y = f(x)
% with `x' in [a, b].
% a - Left integration interval endpoint.
% b - Right integration interval endpoint.
% n - Maximum level in Romberg table.
%
% RETURNS:
% R - Romberg table. Represented as an (n+1)-by-(n+1) lower
% triangular matrix of integral approximations.
%
% SEE ALSO:
% TRAPZ, QUAD, QUADL.

% NOTE: all indices adjusted for MATLAB's one-based indexing scheme.

% Pre-allocate the Romberg table. Avoids subsequent re-allocation which
% is often very costly.
R = zeros([n + 1, n + 1]);

% Initial approximation. Single interval trapezoidal rule.
R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b));

% First column of Romberg table. Increasingly accurate trapezoidal
% approximations.
for i = 1 : n,
h = (b - a) / 2^i;

s = 0;
for k = 1 : 2^(i-1),
s = s + feval(f, a + (2*k - 1)*h);
end

R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s;
end

% Richardson extrapolation gives remainder of Romberg table.
%
% Note: The table is computed by columns rather than the more
% traditional row version. The reason is that this prevents frequent
% (and needless) re-computation of the `fac' quantity.
%
% Moreover, MATLAB matrices internally use ``column major'' ordering so
% this version is less harmful to computer memory cache systems. This
% reason is an implementational detail, though, and less important in
% introctory courses such as MA2501.
for j = 1 : n,
fac = 1 / (4^j - 1);
for m = j : n,
R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1));
end
end

function ff=f(x)
ff=2/sqrt(pi)*exp(-x);
%=================

運行：
>> R=romberg('f', 0, 1, 5)

R =

0.771743332258054 0 0 0 0 0
0.728069946441243 0.713512151168973 0 0 0 0
0.716982762290904 0.713287034240791 0.713272026445579 0 0 0
0.714200167058928 0.713272635314936 0.713271675386546 0.713271669814180 0 0
0.713503839348432 0.713271730111600 0.713271669764711 0.713271669675476 0.713271669674932 0
0.713329714927254 0.713271673453528 0.713271669676323 0.713271669674920 0.713271669674918 0.713271669674918

『貳』 關於matlab，龍貝格法求積分，求高手看一下錯在哪裡

%正解地
function y=Romberg(a,b,e)
n=10;%初始量，大致設一下
h=b-a;
t=zeros(n+1,n+1);
t(1,1)=h/4*(f(a)+2*f((a+b)/2)+f(b));
for i=1:n
t(1,i+1)=1/2*t(1,i);
for x=a+h/(2^i):h/2^(i-1):b-h/(2^i)
t(1,i+1)=t(1,i+1)+h/(2^i)*f(x);
end
end

for i=1:2
t(2,i)=4/3*t(1,i+1)-1/3*t(1,i);
end
i=3;
while(abs(t(2,i-1)-t(2,i-2))>e)
t(2,i)=4/3*t(1,i+1)-1/3*t(1,i);
y=t(2,i);
i=i+1;
if i>n break;
end
end

for i=1:2
t(3,i)=16/15*t(2,i+1)-1/15*t(2,i);
end
i=3;
while(abs(t(3,i-1)-t(3,i-2))>e)
t(3,i)=16/15*t(2,i+1)-1/15*t(2,i);
y=t(3,i);
i=i+1;
if i>n-1 break;
end
end

『叄』 matlab數值計算案例分析的目錄

第1章MATLAB編程基礎1
1.1 矩陣的基本操作與基本運算1
1.1.1 矩陣的基本操作1
1.1.2 矩陣的基本運算2
1.1.3 *與 .*和/與./ 的區別3
1.1.4 使用find函數索引符合某些特定條件的矩陣元素3
1.1.5 eps函數與避免除以0的方法4
1.2 MATLAB的數據結構4
1.3 變數、腳本與函數8
1.3.1 變 量8
1.3.2 全局變數使用例子9
1.3.3 局部變數不會被替代的例子10
1.3.4 函數與腳本10
1.3.5 函數的構成11
1.3.6 函數的類型12
1.3.7 函數調用與函數句柄14
1.3.8 可變參數函數調用14
1.4 MATLAB技巧15
1.4.1 MATLAB的函數重載15
1.4.2 冒號(:)操作符17
1.4.3 Tab鍵自動補全17
1.4.4 上下箭頭回調17
1.4.5 可變參數個數的函數的佔位符17
1.4.6 whos 查看18
1.4.7 whos 通配符的例子18
1.4.8 程序調試18
1.5 MATLAB工具箱函數ode23剖析18
1.6 MATLAB的幫助文檔導航22
1.7 MATLAB常見錯誤23
1.7.1 常見寫法錯誤23
1.7.2 字元串連接出錯24
1.7.3 矩陣維數不同的例子25
1.7.4 賦值出錯26
第2章數值分析的基本概念27
2.1 數值分析的研究對象27
2.2 誤差與有效數字30
2.2.1 誤差的產生及分類30
2.2.2 誤差的相關概念30
2.3 近似計算中的注意事項31
2.4 數值演算法的穩定性34
2.5 機器精度35
第3章數據插值37
3.1 插值與多項式插值37
3.2 Lagrange插值37
3.2.1 Lagrange插值的定義37
3.2.2 Lagrange插值的MATLAB實現38
3.3 Newton插值40
3.3.1 Newton插值定義40
3.3.2 有限差商40
3.3.3 Newton插值的MATLAB實現41
3.4 Hermite插值42
3.4.1 Hermite插值定義42
3.4.2 Hermite插值的MATLAB實現43
3.5 分段低次插值45
3.5.1 高次插值的Runge現象45
3.5.2 分段低次Lagrange插值45
3.5.3 interp1函數46
3.6 三次樣條插值47
3.6.1 三次樣條插值47
3.6.2 三次樣條函數48
第4章數據擬合50
4.1 數據的曲線擬合50
4.1.1 曲線擬合的誤差50
4.1.2 曲線擬合的最小二乘法51
4.2 多項式擬合52
4.2.1 多項式曲線擬合52
4.2.2 多項式曲線擬合的MATLAB實現52
4.2.3 MATLAB多項式曲線擬合應用的擴展54
4.3 圓擬合的例子講解57
4.3.1 圓擬合問題描述(使用最小二乘方法)57
4.3.2 圓擬合的MATLAB實現58
4.4 cftool自定義擬合60
4.5 cftool代碼自動生成與修改62
第5章數值積分66
5.1 數值積分的基本思想66
5.1.1 數值求積的基本思想66
5.1.2 幾種常見的數值積分公式66
5.2 數值求積公式的構造67
5.2.1 代數精度68
5.2.2 插值型求積公式68
5.2.3 Newton-Cotes求積公式69
5.3 復化積分公式70
5.3.1 復化Simpson公式70
5.3.2 復化求積公式及其MATLAB實現70
5.3.3 MATLAB的trapz函數72
5.4 Romberg求積公式73
5.4.1 數值積分公式誤差分析73
5.4.2 Romberg演算法74
5.4.3 Romberg求積公式的MATLAB實現76
5.5 Gauss求積公式77
5.5.1 Gauss積分公式77
5.5.2 Gauss-Legendre求積公式的MATLAB實現及應用實例78
5.6 積分的運算選講79
5.6.1 二重積分79
5.6.2 三重積分79
5.6.3 變上限積分79
5.6.4 符號積分81
5.6.5 MATLAB常見積分函數列表82
第6章常微分方程83
6.1 常微分方程分類及其表示形式83
6.1.1 MATLAB關於ODE的函數幫助簡介83
6.1.2 MATLAB ODE suite中關於ODE的分類83
6.2 典型常微分方程舉例84
6.2.1 一階常微分方程84
6.2.2 二階常微分方程84
6.2.3 高階常微分方程85
6.2.4 邊值問題85
6.2.5 延遲微分方程85
6.3 解的存在性、唯一性和適定性86
6.3.1 初值問題的存在性與唯一性86
6.3.2 MATLAB中常微分方程的通用形式及其向量表示87
6.3.3 剛性常微分方程87
6.4 常微分方程的時域頻域表示以及狀態方程表示89
6.4.1 時域與頻域表示形式89
6.4.2 狀態空間表示形式90
6.5 單步多步和顯式隱式概念91
6.6 常微分方程數值求解方法構造思想舉例92
6.7 常微分方程數值解的基本原理93
6.7.1 一階常微分方程與一階微分方程組93
6.7.2 求解區間[a,b]的離散93
6.7.3 微分方程的離散93
6.7.4 Taylor展開法94
6.7.5 常微分方程數值求解的歐拉方法97
6.7.6 歐拉方法的MATLAB實現97
6.7.7 改進的歐拉方法98
6.7.8 改進的歐拉方法的MATLAB實現99
6.7.9 四階龍格庫塔公式的MATLAB實現99
6.7.10 Adams預測校正公式100
6.8 常微分方程工具箱102
6.8.1 總體介紹102
6.8.2 各個求解器的特點與比較103
6.8.3 使用odefile.m模板求解常微分方程103
6.8.4 odefile.m模板使用105
6.9 單自由度振動系統例子106
6.9.1 單自由度二階系統基於傳遞函數與狀態空間的simulink模型求解106
6.9.2 總 結110
6.10 三自由度振動系統例子110
6.10.1 三自由度振動系統simulink模型求解以及狀態方程的ode45求解器求解110
6.10.2 總 結114
第7章線性方程組的迭代解法115
7.1 線性方程組的迭代法概述115
7.1.1 迭代法概述及壓縮原理115
7.1.2 迭代法基本概念115
7.1.3 MATLAB的相關命令117
7.2 常見的線性方程組的迭代法118
7.2.1 Jacobi迭代法118
7.2.2 Gauss-Seidel迭代法120
7.2.3 SOR迭代法123
7.3 迭代法的收斂性125
7.3.1 迭代法的收斂性定理125
7.3.2 主對角優勢125
7.3.3 SOR迭代法的收斂性126
第7章線性方程組的直接解法127
8.1 線性方程組的消元法127
8.1.1 線性方程組的直接求解方法127
8.1.2 Gauss消去法127
8.1.3 Gauss主元素法130
8.1.4 Jordan消去法133
8.2 矩陣的三角分解135
8.2.1 LU分解136
8.2.2 LU分解的MATLAB實現136
8.2.3 對稱正定矩陣的Cholesky分解138
8.2.4 Cholesky分解法的MATLAB實現139
8.2.5 改進平方根法141
8.2.6 改進平方根法的MATLAB實現142
8.3 MATLAB的相關命令144
8.3.1 逆矩陣144
8.3.2 矩陣的左除及最小二乘解145
8.3.3 欠定方程的解145
第9章非線性方程求解147
9.1 求解非線性方程的MATLAB符號法147
9.2 二分法149
9.2.1 二分法原理149
9.2.2 二分法的MATLAB程序149
9.3 迭代法151
9.3.1 迭代法原理151
9.3.2 迭代法的幾何意義151
9.3.3 迭代法的MATLAB程序152
9.4 切線法153
9.4.1 切線法的幾何意義154
9.4.2 切線法的收斂性154
9.5 割線法（弦截法）155
9.5.1 割線法的幾何意義155
9.5.2 割線法的MATLAB程序155
9.6 常見非線性方程數值方法的優缺點156
9.7 方程f(x)=0數值解的MATLAB實現157
9.7.1 求函數零點指令fzero157
9.7.2 fzero的使用舉例157
9.8 求解非線性方程組MATLAB命令160
9.8.1 符號方程組求解160
9.8.2 求解非線性方程組的基本方法161
9.8.3 求方程組的數值解162
第10章偏微分方程數值解166
10.1 基本概念166
10.2 有限差分法167
10.2.1 橢圓方程的差分形式167
10.2.2 拋物方程的差分形式168
10.2.3 雙曲方程的差分形式170
10.3 MATLAB的pdepe函數171
10.3.1 pdepe函數的說明171
10.3.2 pdepe函數的實例172
10.4 MATLAB的PDEtool工具箱173
10.4.1 PDEtool的界面174
10.4.2 PDEtool的使用174
第11章數值優化177
11.1 單變數函數優化177
11.1.1 基本數學原理177
11.1.2 黃金分割法178
11.1.3 牛頓法181
11.1.4 最速下降法185
11.1.5 共軛梯度法188
11.2 多變數函數優化191
11.2.1 Nelder-mead方法191
11.2.2 Nelder-mead方法的MATLAB實現192
11.2.3 Powell方法193
11.2.4 Powell方法的MATLAB實現194
11.3 MATLAB最優化函數197
11.3.1 MATLAB最優化工具箱介紹197
11.3.2 MATLAB最優化函數介紹198
11.3.3 MATLAB最優化工具介紹201
11.3.4 MATLAB最優化函數應用實例204
第12章特徵值和特徵向量208
12.1 特徵值與特徵向量208
12.1.1 特徵值與特徵向量的定義208
12.1.2 特徵值與特徵向量的計算208
12.1.3 MATLAB的eig命令209
12.2 冪法與反冪法209
12.2.1 冪法的原理210
12.2.2 冪法的MATLAB實現210
12.2.3 反冪法212
12.2.4 反冪法的MATLAB實現213
12.3 對稱矩陣的特徵值——Jacobi方法214
12.3.1 Jacobi方法的原理214
12.3.2 Jacobi方法的MATLAB實現215
12.4 Householder方法217
12.4.1 初等反射矩陣218
12.4.2 用正交相似變換約化矩陣218
12.4.3 演算法的MATLAB實現220
12.5 QR分解與QR方法221
12.5.1 矩陣的QR分解221
12.5.2 計算矩陣特徵值的QR方法222
12.5.3 QR方法的MATLAB實現222
參考文獻224