matlab练习之数值微积分和方程数值求解

Posted on 2020-03-11 02:29  薛向峰  阅读(1450)  评论(0编辑  收藏  举报

 

 

一、符号积分
求符号积分函数:int

格式:int(f,x,a,b)

功能:计算定积分

格式:int(f,x)

功能:计算不定积分

使用int函数之前,先用syms声明x是符号变量

例:

 

代码:

syms x
y1=1/(1+x^4);
fy1=int(y1)
二、数值积分
      在科学研究和工程技术中,经常遇到积分的计算,虽然有些函数的不定积分可以求出其初等函数表示式,但有更多的函数,它们的不定积分不是初等函数,这样就无法利用牛顿莱布尼兹公式求出其定积分,甚至经常遇到只知道函数在一些离散点的值,但函数表达式未知的情况,在上述情况下就必须以数值方法求定积分的近似值。用数值方法求定积分的近似值,通常称为数值积分。

(1) 用梯形法计算积分;
适用于被积函数为离散数据时,求函数的定积分。该函数调用格式:I=trapz(x,y)
例:

clc,clear
format long
ac=@(x)sin(x)./x
x1=pi/4:pi/50:pi/2;
y1=ac(x1);
s1=trapz(x1,y1)
x2=pi/4:pi/100:pi/2;
y2=ac(x2);
s2=trapz(x2,y2)
(2) 基于变步长辛普森法计算积分
该函数调用格式:[I,n]=quad(‘fname’,a,b,Tol,trace)
其中:‘fname是被积函数名
      a,b是积分上下限
      Tol是精度控制值,省却时取0.001
      Trace:控制是否显示展现积分过程,取0不展现
      I:积分值
      n:被积函数调用次数
例:

ac=@(x)sin(x)./x
s=quad(ac,pi/4,pi/2)

(3) 高精度Lobatto积分法,
     格式:z = quadl(Fun,a,b)
(4) 自适应Gauss-Kronrod数值积分
       z = quadgk(Fun,a,b)
(5) 积分法矢量化自适应simpson数值积分
     格式:z = quadv(Fun,a,b)
      一次可以计算多个积分

 

F=@(x,n)1./((1:n)+x.^2);
quadv(@(x)F(x,6),0,1)
三、数值二重积分
 I=dblquad(f,a,b,c,d,tol,trace),求f(x,y)在[a,b]×[c,d]区域上的二重积分。

f=@(x,y)exp(-x.^2/2).*sin(x.^2+y)
I=dblquad(f,-2,2,-1,1)

F=@(x)sqrt(cos(x.^2)+4*sin(4*x.^2)+1);
x1=0:0.01:2*pi;
y1=F(x1);
s1=trapz(x1,y1)
x2=0:0.005:2*pi;
y2=F(x2);
s2=trapz(x2,y2)%梯形法
s3=quad(F,0,pi*2)%基于变步长辛普森法计算积分
---------------------------------
s1 =
   6.7983 + 3.1486i
s2 =
   6.7920 + 3.1557i
s3 =
   6.7992 + 3.1526i

 

F=@(x)log(1+x)./(1+x.^2);
x2=0:0.005:2*pi;
y2=F(x2);
s2=trapz(x2,y2)%梯形法
s3=quad(F,0,pi*2)
------------------------------
s2 =
    0.9995
s3 =
    0.9997

  求解线性方程组

对方程组Ax=b

1 rank(A)=rank(A,b)=n时,无论m=n还是m>n 有唯一解
m=n时,即方程个数等于未知数个数

A=[1 -1 1 -2;2 0 -1 4;3 2 1 0;-1 2 -1 2];
b=[2 4 -1 -4]’;
B=[A,b];
RA=rank(A)
RB=rank(B)
此时RA=RB=4
则X=A\b或X=inv(A)*b(对于本题计算结果是一样的,不过inv要求矩阵为方阵,此外左除的计算速度会快一点)

m>n时,即方程个数大于未知数个数

rank(A)=rank(A,B)<n时 有无穷多解

方法一

(1)、先求通解,null(A)可以解的AX=0的X的值,即X=null(A).
(2)、求解特解,pinv(A)*B=特解。

3 rank(A)不等于rank(A,B)时 无解

此时,x=pinv(A)*B为其最小二乘解

A=[6 5 -2 5;9 -1 4 -1 ;3 4 2 -2 ;3 -9 0 2];
B=[-4;13;1;11];
rank(A)
X=inv(A)*B
---------------------------------
ans =
     4
X =
    0.6667
   -1.0000
    1.5000
   -0.0000

  

Copyright © 2024 薛向峰
Powered by .NET 8.0 on Kubernetes