數值分析matlab上機作業報告

2021-07-25 06:28:40 字數 4107 閱讀 2089

數值求解正方形域上的poisson方程邊值問題

用matlab語言編寫求解此辺值問題的演算法程式,採用下列三種方法,並比較三種方法的計算速度。1、用sor迭代法求解線性方程組au=f,用試演算法確定最佳鬆弛因子;2、用塊 gauss-sediel迭代法求解線性方程組au=f;3、(預條件)共軛斜量法。

用差分代替微分,對poisson方程進行離散化,得到五點格式的線性方程組

寫成矩陣形式au=f。其中

一用sor迭代法求解線性方程組au=f,用試演算法確定最佳鬆弛因子。

1. 基本原理:

gauss-seidel迭代法計算簡單,但是在實際計算中,其迭代矩陣的譜半徑常接近1,因此收斂很慢。為了克服這個缺點,引進乙個加速因子(又稱鬆弛因子)對gauss-seidel方法進行修正加速。

假設已經計算出第k步迭代的解(i=1,2,···,n),要求下一步迭代的解(i=1,2,···,n)。首先,用gauss-seidel迭代格式計算

然後引入鬆弛因子,用鬆弛因子對和作乙個線性組合。

,i=1,2,…,n

將二者合併成為乙個統一的計算公式:

2. 演算法

(1)gauss-seidel迭代法引入鬆弛因子w:

五點格式即為:

(2)計算步驟:

第一步:給鬆弛因子賦初值w=1.1~1.8,給場值u和場源b賦初值

第二步:用不同的w進行迭代計算。置error=0;

計算在計算機上採用動態計算形式

如果|du|>error則error=|du|,如果error第三步:比較不同的w的迭代次數,用kk存放最小迭代次數,用ww和uu存放相應的w及u。

3. 程式

① [u,k]=sor(u,b,w被下面程式呼叫)

%輸入場初值u0、場源b及鬆弛因子w,通過五點差分格式進行迭代運算,

%如果第k+1次的迭代結果與第k次的差小於精度,則可以近似認為第k+1次的迭代

%結果是精確解,然後返回迭代次數k和迭代解

function [u,k]=sor(u,b,w輸出迭代結果u,及迭代次數k

m=length(um為u的維數

h=1/(m-1h為步長

n=10000;e=0.0000001e為精度

for k=1:nk為記錄迭代次數

error=0;

for j=2:m-1

for jj=2:m-1

sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1);

du=w*(h^2*b(jj,j)-sum)/4計算u的修正量

u(jj,j)= u(jj,j)+du修正u

if error end

endif error<=e,break;end判斷是否達到精度

end② [kk,ww,uu]=sor_5dianchafen(n)

%用超鬆弛迭代法求解正方形域上的poisson方程邊值問題

%用5點差分格式求取泊松問題

%輸入n,對x、y軸進行n等分;先確定場u的邊界及場源b,在呼叫[u,k]=sor(u,b,w);

%用不同w計算的迭代次數不同,用kk存放最小的迭代次數,

%用ww和uu分別存放最佳鬆弛因子w和精確解

function [kk,ww,uu]=sor_5dianchafen(n)

w=[1.1:0.1:1.8];m=length(ww為鬆弛因子

kk=1000; ww=0kk是最少迭代次數,ww是最鬆弛因子

h=1/nh步長

b=zeros(n+1,n+1計算場源b

tic;

for i=2:n+1

for j=2:n+1

b(i,j)=(i-1)*(j-1)*h^2;

endenduu=zeros(n+1,n+1); u=zeros(n+1,n+1對u賦初值

u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;

mu=u初值mu以便不同的w計算

for i=1:m用不同的w計算迭代

[u,k]=sor(mu,b,w(i呼叫[u,k]=sor(u,b,w),返回迭代次數及精確解

if kk>k, kk=k;ww=w(i);uu=u;end %把最少迭代次數付給kk,及其w,u賦給ww,uu

endt=toc統計程式運算時間

4.計算結果

>> format short

>> n=10;

>> [kk,ww,uu]=sor_5dianchafen(n)

t = 0.0310

kk =

48ww =

1.6000

uu =

columns 1 through 8

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045

1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088

1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128

1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162

1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189

1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203

1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201

1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175

1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

columns 9 through 11

1.0000 1.0000 1.0000

1.0037 1.0023 1.0000

1.0073 1.0045 1.0000

1.0107 1.0066 1.0000

1.0136 1.0084 1.0000

1.0160 1.0100 1.0000

1.0174 1.0110 1.0000

1.0175 1.0113 1.0000

1.0155 1.0103 1.0000

1.0103 1.0072 1.0000

1.0000 1.0000 1.0000

>> contourf (uu, 'displayname', 'uu'); figure(gcf)

圖一超鬆弛

二用塊jacobi迭代法求解線性方程組au=f。

1. 基本原理:

對a做自然分解a=d-l-u=d-(l+u)

其中d是有a的對角線元素組成的矩陣,l是由a的對角線以下元素組成的矩陣,u是由a得對角線以上元素組成的矩陣。

於是將m=d,n=l+u,代入得到

dx=(l+u)x+b

任取x的初值進行迭代

2. 演算法:

(1)gauss-sediel迭代法原理

五點差分格式:

因為a可以寫成塊狀,即

如果把每一條線上的節點看作乙個組,可以把au=f表示成塊狀求解:

(2)計算步驟:

第一步:給場值u和場源b賦初值,及定義

MATLAB數值計算上機

1 求下列矩陣的主對角元素 上三角陣 下三角陣 逆矩陣 行列式值 秩 範數 條件數 跡 特徵值 特徵向量。2 分別用矩陣求逆 矩陣除法 矩陣的lu分解求下列方程組的解。3 將100個學生5門功課的成績存入矩陣p中,進行如下處理 1 分別求每門課的最高分 最低分及相應學生序號 2 分別求每門課的平均分...

數值分析上機報告

數學與計算科學學院 實驗報告 實驗專案名稱迭代法求線性方程組 所屬課程名稱數值方法a 實驗型別驗證型 實驗日期2011.11.03 班級學號 姓名成績 附錄1 源程式 附錄2 實驗報告填寫說明 1 實驗專案名稱 要求與實驗教學大綱一致。2 實驗目的 目的要明確,要抓住重點,符合實驗教學大綱要求。3 ...

數值分析上機實習報告

學號 姓名 專業 聯絡 2011 年 11 月 序言1.程式語言 我的數值分析上機實踐採用matlab語言作為程式語言,利用matlab集數值計算 符號計算和圖形視覺化三大功能與一體的特點程式設計實現,採用函式檔案編寫程式命令流,並將計算結果儲存為excel格式,並根據計算結果生成圖形,方便查閱 分...