Matlab编程实例

Matlab非线性方程求解器fsolve总结

fsolve是采用最小二乘法来求解非线性方程。它的一般求解方式为：

X=FSOLVE(FUN,X0,OPTIONS)

其中，fun是要求解的非线性方程，X0是变量初值，options由optimset函数产生的结构体，用于对优化参数的设置，可以省略（采用默认值）。

Fsolve可以求解简单的一维非线性方程，如：

x = fsolve(@myfun,[0.5 2 4],optimset('Display','iter')); %求解在初值分别为0.5，2和4时方程的解

其中，函数myfun的定义为：

function F = myfun(x)

F = sin(x);

Fsolve还可以求解大型的非线性方程组，如

x0 = [51.6;rand;unifrnd(-1,1);rand];

h=optimset;

h.MaxFunEvals=20000;

h.MaxIter=5000;

h.Display='off';

[p,fval] = fsolve(@f,x0,options);

此时，方程组可以写成矩阵形式：

function F=f(x)

F=[x(1)+x(2)*(1-exp(-(x(3)*(0)^x(4))))-51.61;

x(1)+x(2)*(1-exp(-(x(3)*(9.78)^x(4))))-51.91;

x(1)+x(2)*(1-exp(-(x(3)*(30.68)^x(4))))-53.27;

x(1)+x(2)*(1-exp(-(x(3)*(59.7)^x(4))))-59.68;];

自编基于龙格库塔法的Matlab数值积分函数

function varargout = rkkt(varargin)

%==============================================

% 采用四阶龙格库塔法数值积分

% rkkt(F,x0,t0,tfinal,ts)

% F为函数名

% x0为积分变量的初始值

% t0为积分初始时刻

% tfinal为积分终止时刻

% ts为积分步长

% example，如果有一个函数名为foo,则，求解指令为

% rkkt(@foo,x0,t0,tfinal,tspan);

% https://www.360docs.net/doc/d210444892.html,/xianfa110

% 可求解任意阶微分方程（组）

% 例如，对于多阶微分方程，我们需对初值进行设定如下：

% rkkt(@foo,[x1_0,x2_0,x3_0,...],t0,tfinal,tspan);

%=================================================== n=length(varargin);

if n==0

disp('请输入函数名！');

return;

end

if n==1

F=varargin{1};

x0=0;

t0=0;

tmax=1;

ts=0.01;

end

if n==2

F=varargin{1};

x0=varargin{2};

t0=0;

tmax=1;

ts=0.01;

end

if n==3

F=varargin{1};

x0=varargin{2};

t0=varargin{3};;

tmax=t0+1;

ts=0.01;

end

if n==4

F=varargin{1};

x0=varargin{2};

t0=varargin{3};

tmax=varargin{4};

ts=0.01;

end

if n==5

F=varargin{1};

x0=varargin{2};

t0=varargin{3};

tmax=varargin{4};

ts=varargin{5};

end

t=t0;

n=length(t:ts:tmax);

m=length(feval_r(F,x0,t0));

x=zeros(n,m);

x(1,:)=x0;

x_temp=x(1,:);

%%%四阶龙格库塔法

for i=2:n

x_0=x_temp;

k1=feval_r(F,x_temp,t);

k1=reshape(k1,1,m);

x_temp=x_0+k1*(ts/2);

t=t+ts/2;

k2=feval_r(F,x_temp,t);

k2=reshape(k2,1,m);

x_temp=x_0+k2*(ts/2);

k3=feval_r(F,x_temp,t);

k3=reshape(k3,1,m);

x_temp=x_0+k3*ts;

t=t+ts/2;

k4=feval_r(F,x_temp,t);

k4=reshape(k4,1,m);

x(i,:)=x_0+ts*(k1+2*k2+2*k3+k4)/6;

x_temp=x(i,:);

end

%%%画图

t=t0:ts:tmax;

length(t);

plot(t,x);grid on;

varargout{1}=x;

自编最小二乘法的Matlab参数辨识程序function [sysd,sys,err] = ID(Y,U,Ts)

%

%基于递推最小二乘法的参数辨识程序

%仅针对二阶系统：)

%出处：https://www.360docs.net/doc/d210444892.html,/xianfa110

%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=

%Inputs:

%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=

%Y = nX1 vector of your model output

%U = nX1 vector of your model input

%Ts = sample time

%

%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=

%Outputs:

%=-=-=-=-=-=-=-=-=-=-==-=-=-=-=-=

%sysd = discrete-time transfer function identified

%sys = continuous-time transfer function identified

%err = error

%

if nargin<3 || nargin>3

error('Must be three inputs!');

end

if length(Y)~=length(U)

error('length of inputs must be equal.');

end

n=length(U);

Y=reshape(Y,n,1);

U=reshape(U,n,1);

theta=[0.1;0.1;0.1;0.1;0.1];

P=2^25*eye(5);

R0=1;

for m=3:n

X=[Y(m-1) Y(m-2) U(m) U(m-1) U(m-2)]';

alfa=1/(R0+X'*P*X);

L=alfa*P*X;

theta(:,m-1)=theta(:,m-2)+L*(Y(m)-X'*theta(:,m-2));

P=P/R0-alfa*P*X*X'*P;

err=Y(m)-X'*theta(:,m-2);

if abs(err)<=1e-10

break;

end

end

m=length(theta(1,:));

result=[-theta(1:2,m);theta(3:5,m)];

t=1:m;

figure;

plot(t,theta(1,:),t,theta(2,:),t,theta(3,:),t,theta(4,:),t,theta(5,:)); legend('th1','th2','th3','th4','th5');

num=[result(3),result(4),result(5)];

den=[1,result(1),result(2)];

sysd=tf(num,den,Ts);

[n,d]=d2cm(num,den,Ts,'tustin');

sys=tf(n,d);

%%====================================================

exaple:

对于以下模型：

运行之后，数据通过Scope传递到工作空间，方法参见Simulink利用Scope输出及绘制仿真波形技巧。

输入以下代码：

Y=data(:,3);U=data(:,2);

[sysd,sys,e]=ID(Y,U,0.001);

得到结果如下：

Transfer function:

-4.314e-009 z^2 + 8.784e-005 z + 4.362e-005

-------------------------------------------

z^2 - 1.975 z + 0.9753

Sampling time: 0.001

Transfer function:

-1.12e-005 s^2 - 0.04417 s + 133.1

----------------------------------

s^2 + 25.02 s - 0.008906

e =

-6.8120e-011

Simulink利用Scope输出及绘制仿真波形技巧

在用Simulink做仿真时，我们经常会用到示波器Scope来观察波形，它可以对波形进行局部放大、按横、纵座标放大，非常方便，但是如果我们要保存波形时，就最好别直接拷贝Scope波形了，因为它的背景是黑的，而且不能进行线形修改和标注，不适合作为文档用图。

一般的做法是将数据输出到工作空间，然后用画图指令Plot画图。输出到工作空间的方法一般有这么几种：

1．添加To Workspace模块；

2．添加out模块；

3．直接用Scope输出。

本人比较懒，一般不再添加其他输出模块，直接选用方法3。当然不是说放一个Scope就能数出数据的，需要对Scope进行设置。设置界面如下：

这里最好把Limit data points to last勾掉，因为很有可能你的数据会超过5000个。勾选Save data to Workspace，变量类型可以选结构体，结构体带时间，以及向量（后面我们会分别介绍这几种变量类型的画图方法）。

运行Simulink，输出完数据，你就可以利用Matlab的画图工具随心所欲的画图了。

下面以一个例子分别介绍三种变量类型的画图方法。

1．输出类型为向量形式。从图上看到，输出了两维时间序列，而实际输出到工作空间的变量ScopeData为三维序列，其中第一列为时间，这正好为我们画图提供了方便。我们可以采用画图命令如下：

figure;

plot(ScopeData(:,1),ScopeData(:,2),'Li