实验7 - 幂法&反幂法&微分方程求解（ode函数）
2018-12-05 本文已影响0人
Xindolia_Ring
- 用反幂法求解下列矩阵的最大特征值以及对应的特征向量，精确到6位数字：
%% invpower.m
%反幂法，精确到6位小数
%A = [6 2 1; 2 3 1; 1 1 1];
function [s,y] = invpower(A, x0, eps, n)
% s为按模最小特征值，y是对应特征向量
k=1;
r=0;
%r相当于λ0
y=x0./max(abs(x0)); %规范化初始向量
[L,U]=lu(A);
z=L\y;
x=U\z;
u=max(x);
s=1/u; % 按模最小为A的逆按模最大的倒数
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);
运行结果
>> 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
- 分别用幂法求下列矩阵的主特征值，反幂法求下列矩阵的模最小特征值，eig函数求全部特征值和特征向量：
%% 幂法
% mipower.m
function [t,y]=mipower(A,x0,eps,n)
% t为所求特征值，y是对应特征向量
k=1;
z=0; %z相当于λ0
y=x0./max(abs(x0)); %规范化初始向量
x=A*y; %迭代格式
b=max(x); %b相当于β
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位小数
%A = [6 2 1; 2 3 1; 1 1 1];
function [s,y] = invpower(A, x0, eps, n)
% s为按模最小特征值，y是对应特征向量
k=1;
r=0;
%r相当于λ0
y=x0./max(abs(x0)); %规范化初始向量
[L,U]=lu(A);
z=L\y;
x=U\z;
u=max(x);
s=1/u; % 按模最小为A的逆按模最大的倒数
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函数
% 全部特征值和特征向量
% E=eig(A)：求矩阵A的全部特征值，构成向量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
运行结果
矩阵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
>>
- 用法介绍
常微分方程数值求解的命令：
求常微分方程的数值解，MATLAB的命令格式为：
[t,y]=solver('odefun',tspan,y0,options)
其中solver选择ode45等函数名，odefun为根据待解方程或方程组编写的m文件名，tspan为自变量的区间[t0,tf]，即准备在那个区间上求解，y0表示初始值，options用于设定误差限制。
命令格式为：
options=odeset('reltol',rt,'abstol',at)
rt输入相对误差，at输入绝对误差。
常用的函数：
| 函数名 | 简介 | 适用对象 |
|---|---|---|
| ode45 | 单步，4/5阶龙格库塔法 | 大部分ODE |
| ode23 | 单步，2/3阶龙格库塔法 | 快速、精度不高的求解 |
| ode113 | 多步，Adams算法 | 误差要求严格或计算复杂 |
- ode23 解非刚性微分方程，低精度，使用Runge-Kutta法的二三阶算法。
- ode45 解非刚性微分方程，中等精度，使用Runge-Kutta法的四五阶算法。
- ode113 解非刚性微分方程，变精度变阶次Adams-Bashforth-Moulton PECE算法。
- 用ode23函数、ode45函数和ode113函数求解下列微分方程并画图比较：【未运行。。。错了不背锅！】
- y'=x+y,y(0)=1, 0<x<3;
- 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
% 调用语句
%[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);
