sor迭代法求matlab程式

1樓：

迭代法 matlab實現**如下

function [x,n] = jacobi(a,b,x0,eps,varargin)

if nargin ==3

eps = 1.0e-6;

m = 200;

elseif nargin<3

disp('輸入引數數目不足3個');

return

elseif nargin ==5

m = varargin;

endd = diag(diag(a)); %%求a的對角矩陣

l = -tril(a,-1); %%求a的下三角矩陣

u = -triu(a,1); %%求a的上三角矩陣

b = d\(l+u);

f = d\b;

x = b*x0+f;

n = 1;%迭代次數

while norm(x-x0)>=eps

x0 = x;

x = b*x0+f

n = n+1;

if(n>=m)

disp('warning:迭代次數太多,可能不收斂!')

return;

endend

執行效果如下：

擴充套件資料：

迭代法的收斂性判別

收斂性判別條件

sor迭代法收斂的充分必要條件是ρ(λω)<1，ρ(λω)與鬆弛因子ω有關。ρ(λω)與ω的關係以及sor方法收斂的條件有如下定理。

定理1：(kahan)對任意的a

，設其對角元皆非零，則對所有實數ω，有：ρ(λω)≥ ω-1。

推論：如果解ax=b的sor方法收斂，則有ω-1<1，即0<ω<2。

定理2：(ostrowski-reich)設a

，a對稱正定，且0<ω<2，則解ax=b的sor方法收斂。

2樓：

1、開啟matlab之後，在命令列視窗中輸入a=[1 2 3 4;5 6 7 8;8 9 2 5;1 2 4 5]，新建一個a方矩陣。

2、在命令列視窗中輸入inv(a)，按回車鍵，可以看到得到了矩陣的逆。

3、使用inv(a)函式求矩陣的逆需要注意的是，a必須是方矩陣，也就是需要行列數相等的矩陣才可以。

4、也可以在命令列視窗輸入help inv，按回車鍵檢視一下inv()函式的用法。

5、在命令列視窗中輸入a^-1，按回車鍵，可以得到矩陣的逆。

6、需要注意的是使用這種方法也需要矩陣為方矩陣或者為標量。

3樓：匿名使用者

function [n,x]=sor22(a,b,x,nm,w,ww)

%用超鬆弛迭代法求解方程組ax=b

%輸入：a為方程組的係數矩陣，b為方程組右端的列向量，x為迭代初值構成的列向量，nm為最大迭代次數，w為誤差精度，ww為鬆弛因子

%輸出：x為求得的方程組的解構成的列向量，n為迭代次數

n=1;

m=length(a);

d=diag(diag(a)); %令a=d-l-u,計算矩陣d

l=tril(-a)+d; %令a=d-l-u,計算矩陣l

u=triu(-a)+d; %令a=d-l-u,計算矩陣u

m=inv(d-ww*l)*((1-ww)*d+ww*u); %計算迭代矩陣

g=ww*inv(d-ww*l)*b; %計算迭代格式中的常數項

%下面是迭代過程

while n<=nm

x=m*x+g; %用迭代格式進行迭代

if norm(x-x,'inf')

disp('迭代次數為');n

disp('方程組的解為');x

return;

%上面：達到精度要求就結束程式，輸出迭代次數和方程組的解

endx=x;n=n+1;

end%下面：如果達到最大迭代次數仍不收斂，輸出警告語句及迭代的最終結果（並不是方程組的解）

disp('在最大迭代次數內不收斂!');

disp('最大迭代次數後的結果為');x

上面是完整的超鬆弛迭代法

如何用matlab編寫的拉格朗日插值演算法的程式、二階龍格-庫塔方法的程式和sor迭代法的程式

4樓：匿名使用者

拉格朗日function y=lagrange(x0,y0,x)

n=length(x0);m=length(x);

for i=1:m

z=x(i);

s=0.0;

for k=1:n

p=1.0;

for j=1:n

if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j));

endend

s=p*y0(k)+s;

endy(i)=s;

end sor迭代法的matlab程式

function [x]=sor_iterative(a,b)

% 用sor迭代求解線性方程組,矩陣a是方陣

x0=zeros(1,length(b)); % 賦初值

tol=10^(-2); % 給定誤差界

n=1000; % 給定最大迭代次數

[n,n]=size(a); % 確定矩陣a的階

w=1; % 給定鬆弛因子

k=1;

% 迭代過程

while k<=n

x(1)=(b(1)-a(1,2:n)*x0(2:n)')/a(1,1);

for i=2:n

x(i)=(1-w)*x0(i)+w*(b(i)-a(i,1:i-1)*x(1:i-1)'-a(i,i+1:n)*x0(i+1:n)')/a(i,i);

endif max(abs(x-x0))<=tol

fid = fopen('sor_iter_result.txt', 'wt');

fprintf(fid,'\n********用sor迭代求解線性方程組的輸出結果********\n\n');

fprintf(fid,'迭代次數: %d次\n\n',k);

fprintf(fid,'x的值\n\n');

fprintf(fid, '%12.8f \n', x);

break;

endk=k+1;

x0=x;

endif k==n+1

fid = fopen('sor_iter_result.txt', 'wt');

fprintf(fid,'\n********用sor迭代求解線性方程組的輸出結果********\n\n');

fprintf(fid,'迭代次數: %d次\n\n',k);

fprintf(fid,'超過最大迭代次數，求解失敗！');

fclose(fid);

end matlab中龍格-庫塔(runge-kutta)方法原理及實現龍格-庫塔(runge-kutta)方法是一種在工程上應用廣泛的高精度單步演算法。由於此演算法精度高，採取措施對誤差進行抑制，所以其實現原理也較複雜。該演算法是構建在數學支援的基礎之上的。

龍格庫塔方法的理論基礎**於泰勒公式和使用斜率近似表達微分，它在積分割槽間多預計算出幾個點的斜率，然後進行加權平均，用做下一點的依據，從而構造出了精度更高的數值積分計算方法。如果預先求兩個點的斜率就是二階龍格庫塔法，如果預先取四個點就是四階龍格庫塔法。一階常微分方程可以寫作：

y'=f(x,y),使用差分概念。

(yn+1-yn)/h= f(xn,yn)推出（近似等於，極限為yn'）

yn+1=yn+h*f(xn,yn)

另外根據微分中值定理，存在0

yn+1=yn+h*f(xn+th,y(xn+th))

這裡k＝f(xn+th,y(xn+th))稱為平均斜率，龍格庫塔方法就是求得k的一種演算法。

利用這樣的原理，經過複雜的數學推導（過於繁瑣省略），可以得出截斷誤差為o(h^5)的四階龍格庫塔公式：

k1＝f(xn,yn);

k2=f(xn+h/2,yn+(h/2)*k1);

k3=f(xn+h/2,yn+(h/2)*k2);

k4=f(xn+h,yn+h*k3);

yn+1=yn+h*(k1+2k2+2k3+k4)*(1/6);

所以，為了更好更準確地把握時間關係，應自己在理解龍格庫塔原理的基礎上，編寫定步長的龍格庫塔函式，經過學習其原理，已經完成了一維的龍格庫塔函式。

仔細思考之後，發現其實如果是需要解多個微分方程組，可以想象成多個微分方程並行進行求解，時間，步長都是共同的，首先把預定的初始值給每個微分方程的第一步，然後每走一步，對多個微分方程共同求解。想通之後發現，整個過程其實很直觀，只是不停的逼近計算罷了。編寫的定步長的龍格庫塔計算函式：

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%參數列順序依次是微分方程組的函式名稱，初始值向量，步長，時間起點，時間終點（引數形式參考了ode45函式）

n=floor((b-a)/h);%求步數

x(1)=a;%時間起點

y(:,1)=y0;%賦初值，可以是向量，但是要注意維數

for ii=1:n

x(ii+1)=x(ii)+h;

k1=ufunc(x(ii),y(:,ii));

k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);

k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);

k4=ufunc(x(ii)+h,y(:,ii)+h*k3);

y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;

%按照龍格庫塔方法進行數值求解

end呼叫的子函式以及其呼叫語句：

function dy=test_fun(x,y)

dy = zeros(3,1);%初始化列向量

dy(1) = y(2) * y(3);

dy(2) = -y(1) + y(3);

dy(3) = -0.51 * y(1) * y(2);

對該微分方程組用ode45和自編的龍格庫塔函式進行比較，呼叫如下：

[t,f] = ode45(@test_fun,[0 15],[1 1 3]);

subplot(121)

plot(t,f)%matlab自帶的ode45函式效果

title('ode45函式效果')

[t1,f1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%測試時改變test_fun的函式維數，別忘記改變初始值的維數

subplot(122)

plot(t1,f1)%自編的龍格庫塔函式效果

title('自編的龍格庫塔函式')

sor迭代法求matlab程式

C程式設計，迭代法解方程，急求用C編寫程式用牛頓迭代法求方程3xxx4xx5x130在x1附近的根，要求精度為10的6次方

高分求一道迭代模型的matlab程式設計問題，題目如圖。答案私信我，千萬私信答得好追加）

求回溯法連續郵資問題的C或C語言程式，急

sor迭代法 求matlab程式

C程式設計，迭代法解方程，急求用C編寫程式用牛頓迭代法求方程3xxx4xx5x130在x1附近的根，要求精度為10的6次方

高分求一道迭代模型的matlab程式設計問題，題目如圖。答案私信我，千萬私信答得好追加）

求回溯法連續郵資問題的C或C語言程式，急

相關推薦

sor迭代法求matlab程式