實(shí)驗(yàn)7 - 冪法&反冪法&微分方程求解（ode函數(shù)）

  1. 用反冪法求解下列矩陣的最大特征值以及對應(yīng)的特征向量，精確到6位數(shù)字：
%% invpower.m
%反冪法，精確到6位小數(shù)
%A = [6 2 1; 2 3 1; 1 1 1];
function [s,y] = invpower(A, x0, eps, n)
% s為按模最小特征值，y是對應(yīng)特征向量
k=1;
r=0;
%r相當(dāng)于λ0
y=x0./max(abs(x0)); %規(guī)范化初始向量
[L,U]=lu(A);
z=L\y;
x=U\z;
u=max(x);
s=1/u; % 按模最小為A的逆按模最大的倒數(shù)
if abs(u-r)<eps %判斷第一次迭代后是否滿足終止條件
    return
end
while abs(u-r)>eps && k<n %終止條件
    k=k+1;
    r=u;
    y=x./max(abs(x));
    z=L\y;
    x=U\z;
    u=max(x);
end
[m,index]=max(abs(x)); %這兩部保證取出的按模最大特征值
s=1/x(index); %是原值而非絕對值
end


%% test1.m
A = [6 2 1; 2 3 1; 1 1 1];
x0 = [1;1;1];
eps = 1e-6;
n=30;
[s,y]=invpower(A, x0, eps, n);

運(yùn)行結(jié)果

>> test1
>> s

s =

    0.5789

>> y

y =

   -0.0461
   -0.3749
    1.0000

>> eig(A)

ans =

    0.5789
    2.1331
    7.2880

>> A*y-s*y

ans =

   1.0e-06 *

   -0.5408
    0.9292
    0.2269
  1. 分別用冪法求下列矩陣的主特征值，反冪法求下列矩陣的模最小特征值，eig函數(shù)求全部特征值和特征向量：
%% 冪法
% mipower.m

function [t,y]=mipower(A,x0,eps,n) 
% t為所求特征值，y是對應(yīng)特征向量
k=1;
z=0; %z相當(dāng)于λ0
y=x0./max(abs(x0)); %規(guī)范化初始向量
x=A*y; %迭代格式
b=max(x); %b相當(dāng)于β
if abs(z-b)<eps %判斷第一次迭代后是否滿足要求
    t=max(x);
    return;
end
while abs(z-b)>eps && k<n
    k=k+1;
    z=b;
    y=x./max(abs(x));
    x=A*y;
    b=max(x);
end
[m,index]=max(abs(x));
t=x(index);
end

%% invpower.m
%反冪法，精確到6位小數(shù)
%A = [6 2 1; 2 3 1; 1 1 1];
function [s,y] = invpower(A, x0, eps, n)
% s為按模最小特征值，y是對應(yīng)特征向量
k=1;
r=0;
%r相當(dāng)于λ0
y=x0./max(abs(x0)); %規(guī)范化初始向量
[L,U]=lu(A);
z=L\y;
x=U\z;
u=max(x);
s=1/u; % 按模最小為A的逆按模最大的倒數(shù)
if abs(u-r)<eps %判斷第一次迭代后是否滿足終止條件
    return
end
while abs(u-r)>eps && k<n %終止條件
    k=k+1;
    r=u;
    y=x./max(abs(x));
    z=L\y;
    x=U\z;
    u=max(x);
end
[m,index]=max(abs(x)); %這兩部保證取出的按模最大特征值
s=1/x(index); %是原值而非絕對值
end

% test2.m
clear;
clc;
%% 初始化
% 矩陣B,C
B = [2 3 2 3;
    3 3 2 -1;
    2 2 4 4 ;
    3 -1 4 4];

C = [4 -1 0 0 0 0;
    -1 4 -1 0 0 0;
    0 -1 4 -1 0 0;
    0 0 -1 4 -1 0;
    0 0 0 -1 4 -1;
    0 0 0 0 -1 4];
% 變量
x1 = [1;1;1;1];
x2 = [1;1;1;1;1;1];
eps = 1e-6;
n=50;

%% 冪法 mipower.m
% 主特征值
[s1,y1]=mipower(B, x1, eps, n);
[s2,y2]=mipower(C, x2, eps, n);
disp('矩陣B主特征值:');
s1
disp('矩陣C主特征值:');
s2
%% 反冪法 invpower.m
% 最小特征值
[s3,y3]=invpower(B, x1, eps, n);
[s4,y4]=invpower(C, x2, eps, n);
disp('矩陣B最小特征值:');
s3
disp('矩陣C最小特征值:');
s4
%% eig函數(shù)
% 全部特征值和特征向量
% E=eig(A)：求矩陣A的全部特征值，構(gòu)成向量E
res1 = eig(B);
res2 = eig(C);
% 直接求矩陣B的特征值和特征向量
[v1,d1] = eig(B,'nobalance');
[v2,d2] = eig(C,'nobalance');
disp('矩陣B的特征值');
v1
disp('矩陣B的特征向量');
d1
disp('矩陣C的特征值');
v2
disp('矩陣C的特征向量');
d2

運(yùn)行結(jié)果

矩陣B主特征值:

s1 =

   10.1930

矩陣C主特征值:

s2 =

   -5.2470

矩陣B最小特征值:

s3 =

   -0.8042

矩陣C最小特征值:

s4 =

    2.1981

矩陣B的特征值

v1 =

   -0.5709    0.6249    0.2618    0.4636
    0.5269   -0.0637    0.7986    0.2838
   -0.3203   -0.7201   -0.0637    0.6122
    0.5421    0.2947   -0.5381    0.5742

矩陣B的特征向量

d1 =

   -2.4951         0         0         0
         0    0.8042         0         0
         0         0    4.4979         0
         0         0         0   10.1930

矩陣C的特征值

v2 =

    0.2319   -0.4179   -0.5211   -0.5211   -0.4179    0.2319
    0.4179   -0.5211   -0.2319    0.2319    0.5211   -0.4179
    0.5211   -0.2319    0.4179    0.4179   -0.2319    0.5211
    0.5211    0.2319    0.4179   -0.4179   -0.2319   -0.5211
    0.4179    0.5211   -0.2319   -0.2319    0.5211    0.4179
    0.2319    0.4179   -0.5211    0.5211   -0.4179   -0.2319

矩陣C的特征向量

d2 =

    2.1981         0         0         0         0         0
         0    2.7530         0         0         0         0
         0         0    3.5550         0         0         0
         0         0         0    4.4450         0         0
         0         0         0         0    5.2470         0
         0         0         0         0         0    5.8019

>>
  • 用法介紹
    常微分方程數(shù)值求解的命令：
    求常微分方程的數(shù)值解，MATLAB的命令格式為：
[t,y]=solver('odefun',tspan,y0,options)

其中solver選擇ode45等函數(shù)名，odefun為根據(jù)待解方程或方程組編寫的m文件名，tspan為自變量的區(qū)間[t0,tf]，即準(zhǔn)備在那個區(qū)間上求解，y0表示初始值，options用于設(shè)定誤差限制。
命令格式為：

options=odeset('reltol',rt,'abstol',at)

rt輸入相對誤差，at輸入絕對誤差。
常用的函數(shù)：

函數(shù)名 簡介 適用對象
ode45 單步，4/5階龍格庫塔法 大部分ODE
ode23 單步，2/3階龍格庫塔法 快速、精度不高的求解
ode113 多步，Adams算法 誤差要求嚴(yán)格或計(jì)算復(fù)雜
  • ode23 解非剛性微分方程，低精度，使用Runge-Kutta法的二三階算法。
  • ode45 解非剛性微分方程，中等精度，使用Runge-Kutta法的四五階算法。
  • ode113 解非剛性微分方程，變精度變階次Adams-Bashforth-Moulton PECE算法。
  1. 用ode23函數(shù)、ode45函數(shù)和ode113函數(shù)求解下列微分方程并畫圖比較：【未運(yùn)行。。。錯了不背鍋！】
  1. y'=x+y,y(0)=1, 0<x<3;
  2. y'=y-2t/y, y(0)=1, 0<t<4;

解：

% dfun1.m
function f1=dfun1(x,y)
f1=x+y;
end

% dfun2.m
function f2=dfun2(t,y)
f2=y-2t./y;
end

% test3.m
%% ode23
[x,y]=ode23(@dfun1,[0 3],1);
[t,y]=ode23(@dfun2,[0 4],1);
plot(x,y);
plot(t,y);

%% ode45
% 調(diào)用語句
%[T,Y] = ode45( @odefun, tspan, y0 );
[x,y]=ode45(@dfun1,[0 3],1);
[t,y]=ode45(@dfun2,[0 4],1);
% 繪圖
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
legend('x','y','z')

%% ode113
[x,y]=ode113(@dfun1,[0 3],1);
[t,y]=ode113(@dfun2,[0 4],1);
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 考試科目：高等數(shù)學(xué)、線性代數(shù)、概率論與數(shù)理統(tǒng)計(jì) 考試形式和試卷結(jié)構(gòu) 一、試卷滿分及考試時(shí)間 試卷滿分為150分，考...
    Saudade_lh閱讀 1,156評論 0贊 0
  • matlab命令 聲明：本文轉(zhuǎn)自：https://www.douban.com/note/136332003/ 侵...
    我就是個初學(xué)者閱讀 14,502評論 0贊 44
  • 哈嘍，大家好 來北京除了玩就是吃吃吃，烤鴨又是排在頭一號，然而動輒排隊(duì)一小時(shí)、兩小時(shí)，甚至還有仨小時(shí)，既心焦又浪費(fèi)...
    齊齊齊齊齊落閱讀 437評論 0贊 0
  • 人切不可操之過急，不然容易出問題。真理之言，有了教訓(xùn)才能更深刻體會此句話之含義。 每每總是在知道自...
    幽幽白書0閱讀 307評論 0贊 1
  • （一） 方曉芬像往常一樣下班回家，洗好澡，貼了張面膜在客廳里看起了連續(xù)劇，她最近在減肥，洗了個蘋果當(dāng)作晚飯。...
    邱有才閱讀 241評論 0贊 0

友情鏈接更多精彩內(nèi)容