feval函数 matlab_基于MATLAB求解三类方程组

基于MATLAB求解三类方程组

第1类 线性方程组

线性方程组是指各个方程变量的最高次幂为1次的方程组。线性方程组主要有定解方程组、不定解方程组、超定方程组和奇异方程组。线性方程组的通解为齐次线性方程组的基础解系加上线性方程组的一个特解。

线性方程组的一般形式为

矩阵形式为  AX=B

其中，A为系数矩阵；B为常数矩阵。

线性方程组解的判定方法

• 如果系数矩阵A的秩不等于增广矩阵[A,B]的秩，那么方程组无解；

• 如果系数矩阵A的秩等于增广矩阵[A,B]的秩且等于方程变量个数n，那么方程组有唯一解；

• 如果系数矩阵A的秩不等于增广矩阵[A,B]的秩且小于方程变量个数n，那么方程组有无穷多解；

采用MATLAB计算的相关基础代码

求矩阵的秩

求秩函数：rank函数

调用格式：

k=rank(A)

k=rank(A,tol)

求方程组特解

x=A\B;

求方程组基础解系

求基础解系函数：null函数

调用格式：

z=null(A)

z=null(A,'r')

1. 定解方程组

定解方程组：方程的个数与方程变量的个数相等，方程组有唯一的解。

示例：求下列定解方程组的解

求解代码：

A=[2,1;2,-3];

B=[1;5];

results=A\B

结果：

x=1,y=-1

此解为定解方程组的唯一解。

2. 不定解方程组

不定解方程组：方程变量的个数多余方程个数的方程组。

示例：求下列不定解方程组的解

求解代码：

A=[2,1,2;2,-3,2];

B=[6;5];

results=A\B

结果：

x=2.875,y=0.25,z=0

此解为不定解方程组的一个特解。

3. 超定方程组

超定方程组：方程变量的个数少于方程个数的方程组。

示例：求下列超定方程组的解

求解代码：

A=[2,1;-2,3;1,4];

B=[3;5;6];

results=A\B

结果：

x=0.2222,y=1.6154

此解为超定方程组的一个最小二乘近似解。

4. 奇异方程组

奇异方程组：方程组系数矩阵为奇异矩阵(行列式等于0)的方程组.

注意：奇异方程组不能直接求解，需增加方程0x+0y=0。

示例：求下列奇异方程组的解

求解代码：

A=[1,-3;2,-6;0,0];

B=[1;-2;0];

results=A\B

结果：

x=1,y=0

此解为奇异方程组的一个解。

第2类 非线性方程组

非线性方程组:指方程组变量中至少有一个的最高次幂高于1次的方程组。

1. 非线性方程的求解

求解函数：fzero函数

调用格式：

x=fzero(fun,x0)

x=fzero(fun,x0,options)

其中，fun为待求解方程; x0为初始值，可为实数标量或2元素实数矢量。options选项属性主要有:Display(显示级别),其值有'off(不显示输出)', 'iter(在每次迭代时显示输出)', 'final(仅显示最终输出)', 'notify(仅在函数为收敛时显示输出)'; FunValCheck(检查目标函数值是否有效),其值有'on(当目标函数返回的值是complex,Inf或NaN时，回显示错误)', 'off(系统默认值)'; OutputFcn; PlotFcns(绘制算法执行过程中的各个进度测量值),其值有@optimplotx(绘制当前点)和@optimplotfval(绘制函数值)；Tolx(关于正标量的终止容差)，默认值为eps(2.2204e-16)

示例：求非线性方程f=sin(cosh(x))的解。

求解代码：

syms x

fun=inline(sin(cosh(x)));

x0=1;

options=optimset('PlotFcns',{@optimplotfval});

x=fzero(fun,x0,options)

box on

set(gca,'fontsize',16,'fontname','times new roman/宋体')

结果：

x=1.8115

2. 非线性方程组的求解

求解函数：fsolve函数

调用格式：

x=fsolve(fun,x0)

x=fsolve(fun,x0,options)

其中，x0为初始值。

示例：求下面非线性方程组的解，其初始值为[0,0]

求解代码：

fun=@root2d;

x0=[0,0];

options=optimoptions('fsolve','Display','none','PlotFcn',@optimplotfirstorderopt);

x=fsolve(fun,x0,options)

box on

set(gca,'fontsize',16,'fontname','times new roman/宋体')

root2d函数

function F=root2d(x)

F(1)=exp(-exp(-(x(1)+x(2))))-x(2)*(1+x(1)^2);

F(2)=x(1)*cos(x(2))+x(2)*sin(x(1))-0.5;

结果：

x1=0.3532,x2=0.6061

第3类 常微分方程组

微分方程是描述动态系统最常用最重要的一类方程。常微分方程是描述一元函数的导数与自变量间关系的方程。

1. 常微分方程组的符号解(解析解)

求解函数：dsolve函数

调用格式：

R=dsolve(‘eq1','eq2',…,‘cond', ‘v’)

其中,eq1,eq2为常微分方程;cond为边界和初始条件;v为指定的变量,若不指定变量,则默认变量变量为t。在表达微分方程时,用字母D表达微分算子,D2和D3分别表示2阶和3阶微分算子,D后所跟的字母为因变量。

示例：求常微分方程组的符号解,其初始值x0=0,y0=0.

求解代码：

syms x y t

[x,y]=dsolve('D2x+2*Dx+x+Dy+y=0','Dx+x+D2y+2*Dy+y=exp(t)','x(0)=0,y(0)=0');

x=simplify(x)

y=simplify(y)

2. 常微分方程组的数值解(近似解)

求解方法：欧拉法(Euler method)和龙格库塔法(Runge-Kutta method)

(1) 欧拉法(Euler method)

欧拉公式：

y_n+1=y_n+hf(x_n,y_n)

改进欧拉公式：

y_f=y_n+hf(x_n,y_n)

y_b=y_n+hf(x_n,y_n)

y_n+1=(y_f+y_b)/2

示例：求常微分方程组的数值解,其初值为x0=0，y0=0。

求解代码：

fcn=inline('y*tan(x)-sec(x)');

x0=0;

y0=0;

xfnl=3;

step=0.05;

[x1,y1]=EulerCAO(fcn,x0,y0,xfnl,step);

[x2,y2]=ModifyEulerCAO(fcn,x0,y0,xfnl,step);

plot(x1,y1,'b-','linewidth',2)

hold on

plot(x2,y2,'r-','linewidth',2)

legend('Euler method','Modify Euler method')

legend boxoff

set(gca,'fontsize',16,'fontname','times new roman')

EulerCAO函数

function [x,y]=EulerCAO(fcn,x0,y0,xfnl,step)

% fcn：常微分方程

% x0：自变量初值

% y0：函数初值

% xfnl：自变量的终值

% step：步长

n=fix((xfnl-x0)/step);

y(1)=y0;

x(1)=x0;

for i=1:n

x(i+1)=x0+i*step;    y(i+1)=y(i)+step*feval(fcn,x(i),y(i));

end

end

ModifyEulerCAO函数

function [x,y]=ModifyEulerCAO(fcn,x0,y0,xfnl,step)

% fcn：常微分方程

% x0：自变量初值

% y0：函数初值

% xfnl：自变量的终值

% step：步长

n=fix((xfnl-x0)/step);

y(1)=y0;

x(1)=x0;

for i=1:n

    x(i+1)=x0+i*step;

    yf=y(i)+step*feval(fcn,x(i),y(i));

    yb=y(i)+step*feval(fcn,x(i+1),yf);

    y(i+1)=(yf+yb)/2;

end

end

结果：

(2)龙格库塔法(Runge-Kutta method)

MATLAB提供的龙格库塔法函数主要有：ode45(单步显式变步长,求解非刚性微分方程),ode23(单步显式变步长,求解非刚性微分方程),ode113(多步，求解非刚性微分方程),ode23t(求解中等刚性微分方程,且无数值阻尼的解),ode15s(多步求解刚性方程显式变步长),ode23s(单步)和ode23tb(隐式，求解刚性微分方程).

示例：求常微分方程组的数值解

求解代码：

a=1;

b=2;

tspan=[0 10];

y0=[0 0.01];

[t,y]=ode45(@(t,y) odefcn(t,y,a,b),tspan,y0);

plot(t,y(:,1),'b-o','linewidth',2)

hold on

plot(t,y(:,2),'r-^','linewidth',2)

set(gca,'fontsize',16,'fontname','times new roman')

odefcn函数

function dydt=odefcn(t,y,a,b)

dydt=zeros(2,1);

dydt(1)=y(2);

dydt(2)=(a/b)*t.*y(1);

结果：

点分享

点收藏

点点赞

点在看

把时间交给阅读