MATLAB

发布日期: 2022-07-04

更新日期: 2022-09-10

文章字数: 2.4k

阅读时长: 11 分

阅读次数:

1.数值微分与数值积分

1).数值微分

a.数值差分和差商的概念

数值差分

数值差商

f(x)在x0点的==微分==接近于函数在该点的差分

f(x)在x0点的==导数==接近于函数在该点的差商

b.向前差分函数

diff(x) 计算向量x的一阶向前差分

diff(x,n) 计算向量x的n阶向前差分

diff(A,n,dim) 计算矩阵A的n阶差分，dim=1：按列计算，dim=2：按行计算

注意向量长度比x少1

举例

x=0:pi/50:2*pi;
y=sin(x);
yy=diff(y)./diff(x);
x=x(1:end-1);
y=cos(x);
plot(x,yy,x,y);
legend('差分求导','cos(x)','Location','NorthEastOutside')

输出示例

差分求导对比

2).数值积分

a.数值积分的实现1----quadl

[l,n]=quad(filename,a,b,tol,trace)

[l,n]=quadl(filename,a,b,tol,trace) 精度高一些

filename: 被积函数名

[a,b]：积分上下限，必须是有限的，不能是无穷大(Inf)

tol：积分精度，默认值为10的-6次幂

trace：是否展现积分过程，非0：展现， 0：不展现，默认为0

返回值l为定积分的值，n为被积函数的调用次数

举例

digits(10);
format long;
f=@(x) 4./(1+x.*x);
[l,n]=quad(f,0,1,10^-8)
[l,n]=quadl(f,0,1,10^-8)

输出示例

l =   3.141592653733437

n =    61

l =   3.141592653589806

n =   48

b.数值积分的实现2----quadgk

[l,err]=quadgk(filename,a,b)

filename：被积函数

[a,b]：积分上下限，可以是无穷大（-Inf,Inf)，也可以是复数

l：定积分值

err：近似误差范围

format long;
f=@(x) 1./(x.*x);
[l,err]=quadgk(f,1,Inf)

输出示例

l =     1

err =     6.314393452555578e-17

c.数值积分的实现----trapz

l=trapz(x,y,dim) 其中x,y是等长的向量

dim=1 按行求积分，dim=2 按列求积分

x=0:1:6;
y=[1,3,4,3,5,3,6];
l=trapz(x,y)
ll=sum(diff(x).*((y(1:end-1)+y(2:end))./2)) %实现方法
plot(x,y,'-rp');
grid on;

输出示例

l =  21.500000000000000

ll =  21.500000000000000

数值积分trapz

d.二重积分的数值解

二重积分

l=integral2(filename,a,b,c,d)

l=quad2d(filename,a,b,c,d)

l=dblquad(filename,a,b,c,d)

f=@(x,y) exp(-x.*x./2).*sin(x.*x+y);
ans1=integral2(f,-2,2,-1,1)
ans2=quad2d(f,-2,2,-1,1)
ans3=dblquad(f,-2,2,-1,1)

输出示例

ans1 =   1.574498159253891

ans2 =   1.574498141468206

ans3 =   1.574493189744943

e.三重积分的数值解

三重积分

l=integral3(filename,a,b,c,d,e,f)

l=triplequad(filename,a,b,c,d,e,f,tol)

f=@(x,y,z) 4.*x.*z.*exp(-z.*z.*y-x.*x);
ans1=integral3(f,0,pi,0,pi,0,1)
ans3=triplequad(f,0,pi,0,pi,0,1,10^-10)

输出示例

ans1 =   1.732762223027684

ans3 =   1.732762222841299

2.线性方程组求解

1).线性方程组的直接解法

a.左除运算符

对于线性方程组Ax=b，则x=A/b

A=[2,1,-5,1;
   1,-5,0,7;
   0,2,1,-1;
   1,6,-1,-4];
b=[13,-9,6,0]';
x=A/b

输出示例

x =
 -66.555555555555543
  25.666666666666664
 -18.777777777777775
  26.555555555555550

b.利用矩阵分解求解

2).线性方程组的迭代解法

a.雅可比迭代法

b.高斯-赛德尔迭代法

3.非线性方程求解

1).fzero函数

x=fzero(filename,x0)

filename：待求根方程左端的函数表达式

x0：初值

执行的是数据搜索过程

f=@(x) x-1./x+5;
ans1=fzero(f,-5)
ans2=fzero(f,1)

x=-10:0.01:10;
y=x-1./x+5;
plot(x,y)
grid on;

输出示例

ans1 =  -5.192582403567252

ans2 =   0.192582403567252

非线性方程求解fzero

2).fsolve函数

x=fsolve(filename,x0,option)

filename：待求根方程左端的函数表达式

x0：初值

option：优化工具箱的优化参数，调用optimset函数完成

f=@(x) x-1./x+5;
ans1=fsolve(f,-5,optimset('Display','Off'))
ans2=fsolve(f,1,optimset('Display','Off'))

x=-10:0.01:10;
y=x-1./x+5;
plot(x,y)
grid on;

输出示例

ans1 =  -5.192582403048145

ans2 =   0.192582403567235

非线性方程求解fzero

fsolve函数可以对方程组进行求解

举例

f=@(x) [sin(x(1))+x(2)+x(3).^2.*exp(x(1)),x(1)+x(2)+x(3),x(1).*x(2).*x(3)];
t=[1,1,1];
ans=fsolve(f,[1,1,1],optimset('Display','Off'))
x=f(ans)

输出示例

ans =   0.022381575402104  -0.022380299910932  -0.000001275491172

x =
   1.0e-06 *
  -0.593079133408355  -0.000000000001941   0.000638901652815

4.函数极值的计算

1).无约束最优化问题

求最小值函数

[xmin,fmin]=fminbnd(filename,x1,x2,option)

[xmin,fmin]=fminsearch(filename,x0,option)

[xmin,fmin]=fminunc(filename,x0,option)

filename：目标函数

x1,x2：左右边界

x0：一个向量，表示极值点的初值

option：优化参数，调用optimset函数

xmin：极小值

fmin：最小值

举例

f=@(x) x-1./x+5;
[xmin,fmin]=fminbnd(f,-10,-1,optimset('Display','Off'))
[xmin,fmin]=fminbnd(f,1,10,optimset('Display','Off'))

输出示例

xmin =  -9.999946678462546

fmin =  -4.899946145244328

xmin =   1.000053455645318

fmin =   5.000106908433283

2).有约束最优化问题

[xmin,fmin]=fmincon(filename,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,option)

filename：目标函数

x0：一个向量，表示极值点的初值

option：优化参数，调用optimset函数

xmin：最小值点

fmin：最小值

其余参数为约束条件，包括线性不等式约束、线性等式约束、x的下界和上界以及定义非线性约束的函数。如果某个约束不存在，则用空矩阵来表示。

举例

有约束规划问题

% fun.m
function f=fun(x)
f=(x(1)-2)^2+(x(2)-1)^2+(x(3)-7)^2+(x(4)-9)^2
end

% yfun.m
function [m,n]=yfun(x)
m=-1*(x(1)^2/4-x(2)^2+x(3)-x(4)^2+1)
n=x(1)^2+x(2)-x(3)+x(4)-99
end

% 脚本文件
x0=[0;0;0;0];
A=[-1,2,-1,1;
   -2,-1,2,-1];
b=[-0.4;-0.5];
lb=[-100;-100;0;0];
ub=[100;100;10;10];
[xmin,fmin]=fmincon(@fun,x0,A,b,[],[],lb,ub,@yfun)

输出示例

xmin =
   9.975955645177839
   0.803897345426697
   7.001160888012390
   5.677572508024824

fmin =  74.692850492068359

5.常微分方程数值求解

1).一般概念

所谓数值解法就是求未知函数在一系列离散点处的近似值

单步法：求y(n+1)时只用到y(n)

多步法：求y(n+1)时，不仅用到y(n)，还要用到y(n-p),其中p=1，2，…，k,k>0的值。

2).求解函数

[t,y]=solver(filename,tspan,y0,option)

solver为使用的函数名。

==注意solver不是一种函数==！！！

filename：函数名

tspan：形式为[t0,tf]，表示求解区间

y0：初始状态向量

option：设置参数

t：时间向量

y：数值解

其中求解函数solver包括：

ode45 非刚性微分方程，中等精度，最主要的求解函数
ode23 非刚性微分方程，低等精度
ode113 非刚性微分方程
ode15s 刚性微分方程和微分代数方程 s:stiff
ode23s 刚性微分方程
ode23t 适度刚醒微分方程和微分代数方程
ode23tb 刚性微分方程
ode15i全隐式微分方程

举例

f=@(t,y) (y^2-t-2)/(4*t+4);
[t,s]=ode23(f,[0,10],[2])

x=0:0.1:10;
y=(x+1).^(1/2)+1;
plot(x,y,'b',t,s,'rp')
legend('标准值','数值解','Location','NorthEastOutside')
axis([0,10,0,10]);

输出示例

常微分方程ode23

使用ode45函数

f=@(t,y) (y^2-t-2)/(4*t+4);
[t,s]=ode45(f,[0,10],[2])

x=0:0.1:10;
y=(x+1).^(1/2)+1;
plot(x,y,'b',t,s,'rp')
legend('标准值','数值解','Location','NorthEastOutside')
axis([0,10,0,10]);

输出示例

常微分方程ode45

该求解函数是为一阶微分方程设计的，求解高阶微分方程，需要先转换为一阶微分方程组。

3).刚性问题

常微分方程解的分量有的变化很快，有的变化很慢，且相差悬殊，这就是所谓的刚性问题。

举例

使用ode45函数

lamda=10^-5;
f=@(t,y) y^2-y^3;
tspan=[0,2/lamda];
tic;
[t,y]=ode45(f,tspan,lamda);
toc;
disp(['计算的点数',num2str(length(t))])

输出示例

历时 0.131789 秒。
计算的点数120565

使用ode23s函数

lamda=10^-5;
f=@(t,y) y^2-y^3;
tspan=[0,2/lamda];
tic;
[t,y]=ode23s(f,tspan,lamda);
toc;
disp(['计算的点数',num2str(length(t))])

输出示例

历时 0.005314 秒。
计算的点数66

6.偏微分方程数值求解

使用pdepe函数

solution=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

solution：x,t相应值对应的偏微分方程数值解

m：平板为0，柱状为1，球面为2

pdefun：定义函数句柄

icfun：初始条件函数句柄

bcfun：边界条件函数句柄

xmesh：向量，用于指定需要针对tspan中每个值求数值解的点

tspan：向量，用于指定需要针对xmesh中每个值求解的点

options：大多数情况下选默认值即可

1).pdefun.m函数文件格式

先将偏微分方程转化为如下形式：

公式1

function [c,f,s]=pdefun(x,t,u,ux)

其中c,f,s均表示上式中的c,f,s，u为上式中的u，ux为u对x的偏导

2).icfun.m函数文件格式

将初值函数转化为如下格式：

公式2

function u0=icfun(x)

其中u0表示初值

3).bcfun.m函数文件格式

将边界条件转化为以下格式：

公式3

function [pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)：

ul是左边界xl=a的近似解，ur是右边界xr=a的近似解。对于左边界pl,ql，需要用ul,对于右边界pr,qr，需要用ur。

4).例题

偏微分例题1

转换形式

偏微分例题分析

% fun.m函数文件
function [c,f,s]=fun(x,t,u,ux)
c=[1;1];
f=[0.024*ux(1);0.017*ux(2)];
y=u(1)-u(2);
F=exp(5.73*y)-exp(-11.46*y);
s=[-F;F];
end

% funstart.m函数文件
function u0=funstart(x)
u0=[1;0];
end

% pdeboundary.m函数文件
function [pl,ql,pr,qr]=pdeboundary(xl,ul,xr,ur,t)
pl=[0;ul(2)];
ql=[1;0];
pr=[ur(1)-1;0];
qr=[0;1];
end

%脚本文件
x=linspace(0,5,100);
t=linspace(10,20,100);
s=pdepe(0,@fun,@funstart,@pdeboundary,x,t);

figure;
surf(x,t,s(:,:,1));
xlabel('x');
ylabel('t');
zlabel('s');
title('u1');
figure;
surf(x,t,s(:,:,2));
xlabel('x');
ylabel('t');
zlabel('s');
title('u2');

输出示例

偏微分例题u1

偏微分例题u2

易安

https://jingblog.tech/2022/07/04/shu-zhi-wei-ji-fen-yu-fang-cheng-qiu-jie/

本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源易安 !

MATLAB 数学计算

ubuntu20.04防火墙设置

2022-08-22 网站设计

网站服务器防火墙

数据分析与多项式计算

2022-04-03 MATLAB

MATLAB 数学计算

数值微积分与方程求解

1.数值微分与数值积分

1).数值微分

a.数值差分和差商的概念

b.向前差分函数

2).数值积分

a.数值积分的实现1----quadl

b.数值积分的实现2----quadgk

c.数值积分的实现----trapz

d.二重积分的数值解

e.三重积分的数值解

2.线性方程组求解

1).线性方程组的直接解法

a.左除运算符

b.利用矩阵分解求解

2).线性方程组的迭代解法

a.雅可比迭代法

b.高斯-赛德尔迭代法

3.非线性方程求解

1).fzero函数

2).fsolve函数

4.函数极值的计算

1).无约束最优化问题

2).有约束最优化问题

5.常微分方程数值求解

1).一般概念

2).求解函数

3).刚性问题

6.偏微分方程数值求解

1).pdefun.m函数文件格式

2).icfun.m函数文件格式

3).bcfun.m函数文件格式

4).例题