matlab动力学分析程序详解

1

1.微分方程的定义

对于duffing 方程03

2

=++x x x ω

，先将方程写作???

--==3

1122

21x x x x x ω function dy=duffing(t,x) omega=1;%定义参数 f1=x(2);

f2=-omega^2*x(1)-x(1)^3; dy=[f1;f2];

2.微分方程的求解

function solve (tstop) tstop=500;%定义时间长度 y0=[0.01;0];%定义初始条件

[t,y]=ode45('duffing',tstop,y0,[]);

function solve (tstop) step=0.01;%定义步长

y0=rand(1,2);%随机初始条件

tspan=[0:step:500];%定义时间范围 [t,y]=ode45('duffing',tspan,y0);

3.时间历程的绘制

时间历程横轴为t ，纵轴为y ，绘制时只取稳态部分。 plot(t,y(:,1));%绘制y 的时间历程 xlabel('t')%横轴为t ylabel('y')%纵轴为y grid;%显示网格线

2

axis([460 500 -Inf Inf])%图形显示范围设置

4.相图的绘制

相图的横轴为y ，纵轴为dy/dt ，绘制时也只取稳态部分。红色部分表示只取最后1000个点。

plot(y(end-1000:end ,1),y(end-1000:end ,2));%绘制y 的时间历程

xlabel('y')%横轴为y

ylabel('dy/dt')%纵轴为dy/dt grid;%显示网格线

5.Poincare 映射的绘制

对于不同的系统，Poincare 截面的选取方法也不同

对于自治系统一般每过其对应线性系统的固有周期，截取一次 对于非自治系统，一般每过其激励的周期，截取一次

例程：duffing 方程03

2=++x x x ω

的poincare 映射 function poincare(tstop)

global omega; omega=1;

T=2*pi/omega;%线性系统的周期或激励的周期 step=T/100;%定义步长为T/100 y0=[0.01;0];%初始条件

tspan=[0:step:100*T];%定义时间范围 [t,y]=ode45('duffing',tspan,y0);

for i=5000:100:10000%稳态过程每个周期取一个点 plot(y(i,1),y(i,2),'b.'); hold on;% 保留上一次的图形 end

xlabel('y');ylabel('dy/dt');

3

Poincare 映射也可以通过取极值点得到 function poincare(tstop) y0=[0.01;0];

tspan=[0:0.01:500];

[t,y]=ode45('duffing',tspan,y0); count=find(t>100);%截取稳态过程 y=y(count,:);

n=length(y(:,1));%计算点的总数 for i=2:n-1

if y(i-1,1)+epsy(i+1,1)+eps % 简单的取出局部最大值

plot(y(i,1),y(i,2),'.'); hold on end end

xlabel('y');ylabel('dy/dt');

6.频谱

yy=fft(y(end-1000:end,1)); N=length(yy); power=abs(yy);

freq=(1:N-1)*1/step/N;

plot(freq(1:N/2),power(1:N/2)); xlabel('f(y)') ylabel('y')

7.算例

duffing 方程03

=++x x x

的时间历程，相图，频谱和poincare 映射。

function dy=duffing(t,x)

omega=1;%定义参数

f1=x(2);

f2=-omega^2*x(1)-x(1)^3;

dy=[f1;f2]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function duffsim(tstop)

step=0.01

y0=[0.1;0];

tspan=[0:step:500];

[t,y]=ode45('duffing',tspan,y0); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,1)

plot(t,y(:,1));%绘制y的时间历程

xlabel('t')%横轴为t

ylabel('y')%纵轴为y

grid;%显示网格线

axis([460 500 -Inf Inf])%显示范围设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,2)

plot(y(end-1000:end,1),y(end-1000:end,2));%绘制y的时间历程

xlabel('y')%横轴为y

ylabel('dy/dt')%纵轴为dy/dt

grid;%显示网格线%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,3)

yy=fft(y(end-1000:end,1));

4

N=length(yy);

power=abs(yy);

freq=(1:N-1)*1/step/N;

plot(freq(1:N/2),power(1:N/2));

xlabel('f(y)')

ylabel('y') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,4)

count=find(t>100);%截取稳态过程

y=y(count,:);

n=length(y(:,1));%计算点的总数

for i=2:n-1

if y(i-1,1)+epsy(i+1,1)+eps % 简单的取出局部最大值

plot(y(i,1),y(i,2),'.');hold on;

end

end

xlabel('y');ylabel('dy/dt');

5

8.分岔图的绘制

随F变化的分岔图。

+

+

3.03=

-

x2.1

F

t

x

x

x

cos

function dy=duffing(t,x)

global c;

omega=1;%定义参数

f1=x(2);

f2=omega^2*x(1)-x(1)^3-0.3*x(2)+c*cos(1.2*t); dy=[f1;f2]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

global c; %定义全局变量

range=[0.1:0.002:0.9];%定义参数变化范围

k=0;

6

YY=[];%定义空数组

for c=range

y0=[0.1;0];%初始条件

k=k+1;

tspan=[0:0.01:400];

[t,Y]=ode45('duffing',tspan,y0);

count=find(t>200);

Y=Y(count,:);

j=1;

n=length(Y(:,1));

for i=2:n-1

if Y(i-1,1)+epsY(i+1,1)+eps % 简单的取出局部最大值。

YY(k,j)=Y(i,1);

j=j+1;

end

end

if j>1

plot(c,YY(k,[1:j-1]),'k.','markersize',3);

end

hold on;

index(k)=j-1;

end

xlabel('c');

ylabel('y');

7

随F变化的分岔图

F=0.20

8

9

F=0.27

F=0.275

F=0.2875

10

F=0.32

F=0.36

F=0.4

11

F=0.652

F=0.8

现代信号处理Matlab仿真——例611

例6.11 利用卡尔曼滤波估计一个未知常数 题目： 设已知一个未知常数x 的噪声观测集合，已知噪声v(n)的均值为零， 方差为 ，v(n)与x 不相关，试用卡尔曼滤波估计该常数 题目分析： 回忆Kalman 递推估计公式 由于已知x 为一常数，即不随时间n 变化，因此可以得到： 状态方程： x(n)=x(n-1) 观测方程： y(n)=x(n)+v(n) 得到A(n)=1,C(n)=1, , 将A(n)=1,代入迭代公式 得到：P(n|n-1)=P(n-1|n-1) 用P(n-1)来表示P(n|n-1)和P(n-1|n-1),这是卡尔曼增益表达式变为 从而 2v σ1??(|1)(1)(1|1)(|1)(1)(1|1)(1)()()(|1)()[()(|1)()()]???(|)(|1)()[()()(|1)](|)[()()](|1)H w H H v x n n A n x n n P n n A n P n n A n Q n K n P n n C n C n P n n C n Q n x n n x n n K n y n C n x n n P n n I K n C n P n n --=----=----+=--+=-+--=--2()v v Q n σ=()0w Q n =(|1)(1)(1|1)(1)()H w P n n A n P n n A n Q n -=----+21 ()(|1)[(|1)]v K n P n n P n n σ-=--+22(1)()[1()](1)(1)v v P n P n K n P n P n σσ-=--=-+

matlab在统计数据的描述性分析的应用

统计数据的描述性分析 一、实验目的 熟悉在matlab中实现数据的统计描述方法，掌握基本统计命令：样本均值、样本中位数、样本标准差、样本方差、概率密度函数pdf、概率分布函数df、随机数生成rnd。 二、实验内容 1 、频数表和直方图 数据输入,将你班的任意科目考试成绩输入 >> data=[91 78 90 88 76 81 77 74]; >> [N,X]=hist(data,5) N = 3 1 1 0 3 X = 75.7000 79.1000 82.5000 85.9000 89.3000 >> hist(data,5)

2、基本统计量 1) 样本均值 语法: m=mean(x) 若x 为向量，返回结果m是x 中元素的均值； 若x 为矩阵，返回结果m是行向量，它包含x 每列数据的均值。 2) 样本中位数 语法: m=median(x) 若x 为向量，返回结果m是x 中元素的中位数； 若x 为矩阵，返回结果m是行向量，它包含x 每列数据的中位数3) 样本标准差 语法:y=std(x) 若x 为向量，返回结果y 是x 中元素的标准差； 若x 为矩阵，返回结果y 是行向量，它包含x 每列数据的标准差

std(x)运用n-1 进行标准化处理，n是样本的个数。 4) 样本方差 语法：y=var(x); y=var(x,1) 若x 为向量，返回结果y 是x 中元素的方差； 若x 为矩阵，返回结果y 是行向量，它包含x 每列数据的方差 var(x)运用n-1 进行标准化处理（满足无偏估计的要求），n 是样本的个数。var(x,1)运用n 进行标准化处理，生成关于样本均值的二阶矩。 5) 样本的极差（最大之和最小值之差） 语法：z= range(x) 返回结果z是数组x 的极差。 6) 样本的偏度 语法：s=skewness(x) 说明：偏度反映分布的对称性，s>0 称为右偏态，此时数据位于均值右边的比左边的多；s<0,情况相反；s 接近0 则可认为分布是对称的。 7) 样本的峰度 语法：k= kurtosis(x) 说明：正态分布峰度是3，若k 比3 大得多，表示分布有沉重的尾巴，即样本中含有较多远离均值的数据，峰度可以作衡量偏离正态分布的尺度之一。 >> mean(data) ,

MATLAB数据分析与多项式计算(M)

第7章 MATLAB数据分析与多项式计算 6.1 数据统计处理 6.2 数据插值 6.3 曲线拟合 6.4 离散傅立叶变换 6.5 多项式计算 6.1 数据统计处理 6.1.1 最大值和最小值 MATLAB提供的求数据序列的最大值和最小值的函数分别为max 和min，两个函数的调用格式和操作过程类似。 1．求向量的最大值和最小值 求一个向量X的最大值的函数有两种调用格式，分别是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。 (2) [y,I]=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。 求向量X的最小值的函数是min(X)，用法和max(X)完全相同。 例6-1 求向量x的最大值。 命令如下： x=[-43,72,9,16,23,47]; y=max(x) %求向量x中的最大值 [y,l]=max(x) %求向量x中的最大值及其该元素的位置 2．求矩阵的最大值和最小值 求矩阵A的最大值的函数有3种调用格式，分别是： (1) max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i 列上的最大值。 (2) [Y,U]=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。 (3) max(A,[],dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。 求最小值的函数是min，其用法和max完全相同。

例6-2 分别求3×4矩阵x中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。 3．两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： (1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B 同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 (2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 min函数的用法和max完全相同。 例6-3 求两个2×3矩阵x, y所有同一位置上的较大元素构成的新矩阵p。 6.1.2 求和与求积 数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为： sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘积。 sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。 prod(A)：返回一个行向量，其第i个元素是A的第i列的元素乘积。 sum(A,dim)：当dim为1时，该函数等同于sum(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。 prod(A,dim)：当dim为1时，该函数等同于prod(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。 例6-4 求矩阵A的每行元素的乘积和全部元素的乘积。 6.1.3 平均值和中值 求数据序列平均值的函数是mean，求数据序列中值的函数是median。两个函数的调用格式为： mean(X)：返回向量X的算术平均值。 median(X)：返回向量X的中值。

Matlab仿真实例-卫星轨迹

卫星轨迹 一．问题提出 设卫星在空中运行的运动方程为： 其中是k 重力系数(k=401408km3/s)。卫星轨道采用极坐标表示，通过仿真，研究发射速度对卫星轨道的影响。实验将作出卫星在地球表面（r=6400KM ，θ=0）分别以v=8KM/s,v=10KM/s,v=12KM/s 发射时，卫星绕地球运行的轨迹。 二．问题分析 1．卫星运动方程一个二阶微分方程组，应用Matlab 的常微分方程求解命令ode45求解时，首先需要将二阶微分方程组转换成一阶微分方程组。若设，则有： 2．建立极坐标如上图所示，初值分别为：卫星径向初始位置，即地球半径：y(1,1)=6400；卫星初始角度位置：y(2,1)=0；卫星初始径向线速度：y(3,1)=0；卫星初始周向角速度：y(4,1)=v/6400。 3．将上述一阶微分方程及其初值带入常微分方程求解命令ode45求解，可得到一定时间间隔的卫星的径向坐标值y(1)向量；周向角度坐标值y(2)向量；径向线速度y(3)向量；周向角速度y(4)向量。 4．通过以上步骤所求得的是极坐标下的解，若需要在直角坐标系下绘制卫星的运动轨迹，还需要进行坐标变换，将径向坐标值y(1)向量；周向角度坐标值y(2)向量通过以下方程转换为直角坐标下的横纵坐标值X,Y 。 5．卫星发射速度速度的不同将导致卫星的运动轨迹不同，实验将绘制卫星分别以v=8KM/s ，v=10KM/s ，v=12KM/s 的初速度发射的运动轨迹。 三．Matlab 程序及注释 1．主程序 v=input('请输入卫星发射速度单位Km/s ：\nv=');%卫星发射速度输入。 axis([-264007000-1000042400]);%定制图形输出坐标范围。 %为了直观表达卫星轨迹，以下语句将绘制三维地球。 [x1,y1,z1]=sphere(15);%绘制单位球。 x1=x1*6400;y1=y1*6400;???????-=+-=dt d dt dr r dt d dt d r r k dt r d θ θθ2)(2 22222θ==)2(,)1(y r y ?????????????**-=**+*-===)1(/)4()3(2)4()4()4()1()1()1()3()4()2() 3()1(y y y dt dy y y y y y k dt dy y dt dy y dt dy ???*=*=)] 2(sin[)1(Y )]2(cos[)1(X y y y y

Matlab对采样数据进行频谱分析

使用Matlab对采样数据进行频谱分析 1、采样数据导入Matlab 采样数据的导入至少有三种方法。 第一就是手动将数据整理成Matlab支持的格式，这种方法仅适用于数据量比较小的采样。 第二种方法是使用Matlab的可视化交互操作，具体操作步骤为：File --> Import Data，然后在弹出的对话框中找到保存采样数据的文件，根据提示一步一步即可将数据导入。这种方法适合于数据量较大，但又不是太大的数据。据本人经验，当数据大于15万对之后，读入速度就会显著变慢，出现假死而失败。 第三种方法，使用文件读入命令。数据文件读入命令有textread、fscanf、load 等，如果采样数据保存在txt文件中，则推荐使用 textread命令。如 [a,b]=textread('data.txt','%f%*f%f'); 这条命令将data.txt中保存的数据三个三个分组，将每组的第一个数据送给列向量a，第三个数送给列向量b，第二个数据丢弃。命令类似于C语言，详细可查看其帮助文件。文件读入命令录入采样数据可以处理任意大小的数据量，且录入速度相当快，一百多万的数据不到20秒即可录入。强烈推荐！ 2、对采样数据进行频谱分析 频谱分析自然要使用快速傅里叶变换FFT了，对应的命令即 fft ，简单使用方法为：Y=fft(b,N)，其中b即是采样数据，N为fft数据采样个数。一般不指定N，即简化为Y=fft(b)。Y即为FFT变换后得到的结果，与b的元素数相等，为复数。以频率为横坐标，Y数组每个元素的幅值为纵坐标，画图即得数据b的幅频特性；以频率为横坐标，Y数组每个元素的角度为纵坐标，画图即得数据b的相频特性。典型频谱分析M程序举例如下： clc fs=100; t=[0:1/fs:100]; N=length(t)-1;%减1使N为偶数 %频率分辨率F=1/t=fs/N p=1.3*sin(0.48*2*pi*t)+2.1*sin(0.52*2*pi*t)+1.1*sin(0.53*2*pi*t)... +0.5*sin(1.8*2*pi*t)+0.9*sin(2.2*2*pi*t); %上面模拟对信号进行采样，得到采样数据p，下面对p进行频谱分析 figure(1) plot(t,p); grid on title('信号 p(t)'); xlabel('t') ylabel('p')

数值分析的matlab实现

第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *

数值分析的MATLAB程序

列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);

end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')

MATLAB实现通信系统仿真实例

补充内容：模拟调制系统的MATLAB 仿真 1.抽样定理 为了用实验的手段对连续信号分析，需要先对信号进行抽样（时间上的离散化），把连续数据转变为离散数据分析。抽样（时间离散化）是模拟信号数字化的第一步。 Nyquist 抽样定律：要无失真地恢复出抽样前的信号，要求抽样频率要大于等于两倍基带信号带宽。 抽样定理建立了模拟信号和离散信号之间的关系，在Matlab 中对模拟信号的实验仿真都是通过先抽样，转变成离散信号，然后用该离散信号近似替代原来的模拟信号进行分析的。 【例1】用图形表示DSB 调制波形)4cos()2cos(t t y ππ= 及其包络线。 clf %%计算抽样时间间隔 fh=1;%%调制信号带宽(Hz) fs=100*fh;%%一般选取的抽样频率要远大于基带信号频率，即抽样时间间隔要尽可能短。 ts=1/fs; %%根据抽样时间间隔进行抽样，并计算出信号和包络 t=(0:ts:pi/2)';%抽样时间间隔要足够小，要满足抽样定理。 envelop=cos(2*pi*t);%%DSB 信号包络 y=cos(2*pi*t).*cos(4*pi*t);%已调信号 %画出已调信号包络线 plot(t,envelop,'r:','LineWidth',3); hold on plot(t,-envelop,'r:','LineWidth',3); %画出已调信号波形 plot(t,y,'b','LineWidth',3); axis([0,pi/2,-1,1])% hold off% xlabel('t'); %写出图例 【例2】用图形表示DSB 调制波形)6cos()2cos(t t y ππ= 及其包络线。 clf %%计算抽样时间间隔 fh=1;%%调制信号带宽(Hz) fs=100*fh;%抽样时间间隔要足够小，要满足抽样定理。 ts=1/fs; %%根据抽样时间间隔进行抽样

Matlab大数据处理

Matlab大数据处理2：硬盘访问.mat文件 分类：Matlab Hack2013-09-08 20:16 146人阅读评论(0) 收藏举报Matlab程序中经常要访问.mat文件，通常在作法是用load函数直接加载.mat文件。如果.mat文件非常大，超过了系统可用内存的时候该怎么办呢？Matlab2013b为提供了matfile函数，matfile函数可以通过索引直接访问.mat文件中的Matlab变量，而无需将.mat文件加载入内存。 matfile有两种用法： m = matfile(filename)，用文件名创建matfile对象，通过这个对象可以直接访问mat文件中的matlab变量。 m = matfile(filename,'Writable',isWritable)，isWritable开启或关闭文件写操作。 使用示例： 1. 向mat文件中写入变量 x = magic(20); m = matfile('myFile.mat'); % 创建一个指向myFile.mat的matfile对象 m.x = x; % 写入x m.y(81:100,81:100) = magic(20); % 使用坐标索引

2. 加载变量 filename = 'topography.mat'; m = matfile(filename); topo = m.topo; %读取变量topo [nrows,ncols] = size(m,'stocks'); %读取stocks变量的size avgs = zeros(1,ncols); for idx = 1:ncols avgs(idx) = mean(m.stocks(:,idx)); end 3. 开启写权限 filename = 'myFile.mat'; m = matfile(filename,'Writable',true); 或者 m.Properties.Writable = true;

数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序

（Hilbert 矩阵）病态线性方程组的求解 理论分析表明，数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ?=，11 ij h i j = +-，i ，j = 1,2，…，n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n ，取(1,1,,1)n x =∈K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭 代，SOR 迭代和共轭梯度法求解，比较结果。 3. 结合计算结果，试讨论病态线性方程组的求解。 第1小题： condition.m %第1小题程序 t1=20;%阶数n=20 x1=1:t1; y1=1:t1; for i=1:t1 H=hilb(i); y1(i)=log(cond(H)); end plot(x1,y1); xlabel('阶数n'); ylabel('2-条件数的对数（log(cond(H)）'); title('2-条件数的对数（log(cond(H)）与阶数n 的关系图'); t2=200;%阶数n=200 x2=1:t2; y2=1:t2; for i=1:t2 H=hilb(i); y2(i)=log(cond(H)); end plot(x2,y2); xlabel('阶数n'); ylabel('2-条件数的对数（log(cond(H)）'); title('2-条件数的对数（log(cond(H)）与阶数n 的关系图'); 画出Hilbert 矩阵2-条件数的对数和阶数的关系

n=200时 n=20时 从图中可以看出， 1）在n小于等于13之前，图像近似直线 log（cond（H））~1.519n-1.833 2）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势 第2小题： solve.m%m第2小题主程序 N=4000;

同济大学数值分析matlab编程题汇编

MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8

数值分析算法在matlab中的实现

数值分析matlab实现高斯消元法： function[RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意：因为RA=RB0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1

第6章matlab数据分析与多项式计算_习题答案

第6章 MATLAB数据分析与多项式计算 习题6 一、选择题 1．设A=[1,2,3,4,5;3,4,5,6,7]，则min(max(A))的值是（）。B A．1 B．3 C．5 D．7 2．已知a为3×3矩阵，则运行mean(a)命令是（）。B A．计算a每行的平均值 B．计算a每列的平均值 C．a增加一行平均值 D．a增加一列平均值 3．在MATLAB命令行窗口输入下列命令： >> x=[1,2,3,4]; >> y=polyval(x,1); 则y的值为（）。 D A．5 B．8 C．24 D．10 4．设P是多项式系数向量，A为方阵，则函数polyval(P,A)与函数polyvalm(P,A)的值（）。D A．一个是标量，一个是方阵 B．都是标量 C．值相等 D．值不相等 5．在MATLAB命令行窗口输入下列命令： >> A=[1,0,-2]; >> x=roots(A); 则x(1)的值为（）。 C A．1 B．-2 C． D． 6．关于数据插值与曲线拟合，下列说法不正确的是（）。A A．3次样条方法的插值结果肯定比线性插值方法精度高。 B．插值函数是必须满足原始数据点坐标，而拟合函数则是整体最接近原始数据点，而不一定要必须经过原始数据点。 C．曲线拟合常常采用最小二乘原理，即要求拟合函数与原始数据的均方误差达到极小。 D．插值和拟合都是通过已知数据集来求取未知点的函数值。 二、填空题 1．设A=[1,2,3;10 20 30;4 5 6]，则sum(A)= ，median(A)= 。 [15 27 39]，[4 5 6[ 2．向量[2,0,-1]所代表的多项式是。2x2-1 3．为了求ax2+bx+c=0的根，相应的命令是（假定a、b、c已经赋值）。为了

数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法')

a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2)0 a=p; else b=p; end end if i==n0&&~(fp==0||(abs((b-a)/2)

实验一数据处理方法MATLAB实现

实验一数据处理方法的MATLAB实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容 1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据：实验一数据处理方法的MATLAB实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容

1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据： x=[1,4,2,81,23,45] x = 1 4 2 81 2 3 45 单击保存按钮，保存在Matlab指定目录（C:\Program Files\MATLAB71）下，文件名为“mydat.mat”。 2. 在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 继续在命令窗口中输入命令： （1）求取最大值“max（a）”； >> max(x) ans = 81 （2）求取最小值“min（a）”； >> min(x) ans = 1 （3）求取均值“mean（a）”； >> mean(x) ans =

数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286

B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1)

数值分析 matlab 实验4

(1) 解题过程如下： （1）MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件：1）复化梯形公式文件： function s=T_fuhua(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; 2）复化辛普森公式文件： function s=S_fuhua(f,a,b,n) h=0; h=(b-a)./(2*n); s1=0; https://www.wendangku.net/doc/5216590029.html, -5- s2=0; for k=1:n-1 x=a+h*2*k; s1=s1+feval(f,x); end for k=1:n x=a+h*(2*k-1); s2=s2+feval(f,x); end

s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3; 在MATLAB中输入： f=inline('x/(4+x^2)');a=0;b=1; %inline 构造内联函数对象 for n=2:10 s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10); %调用复化梯形公式，生成任意精度的数值 end exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10) %求出积分的精确值 输出结果：exact = .1115717755 s = Columns 1 through 6 0.1088 0.1104 0.1109 0.1111 0.1113 0.1114 Columns 7 through 9 0.1114 0.1114 0.1115 在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线： r=abs(exact-s); n=2:10; plot(double(n),double(r(n-1))) 得到结果如图所示： （2）在 MATLAB中输入以下程序代码： f=inline('x/(4+x^2)');a=0;b=1;n=9; %inline 构造内联函数对象 t=T_fuhua(f,a,b,n);t=vpa(t,10) s=S_fuhua(f,a,b,n);s=vpa(s,10)

数值分析实验— MATLAB实现

数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月

一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13)

一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5

数值分析五个题目的C语言及Matlab程序

本文档包含上一个文档中的五个数值分析实验题C语言程序及Matlab程序实验一 C程序 #include "stdio.h" #include "math.h" void main() { inti=0; float a=0.1,b=1.9,t=0.0,e=1.9; if((pow(a,7)-28*pow(a,4)+14)*(pow(b,7)-28*pow(b,4)+14)<0) if((7*pow(x,6)-112*pow(x,3))) printf("x=%f,i=%d,e=%f\n",x,i,e); for(i=1;i<7&&e>0.00001;i++) { t=x; x=x-(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3)); e=fabs(t-x); printf("x=%f,i=%d,e=%f\n",x,i,e);

} } Matable 程序 i=0; x=1.9;t=0.0;e=1.9; disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); for i=1:7 t=x; x=x-(x^7-28*x^4+14)/(7*x^6-112*x^3); e=abs(t-x); disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); if e<0.00001 break; end end 实验二 C程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10};