matlab数值微积分与方程数值求解

电子一班王申江

实验九数值微积分与方程数值求解

一、实验目的

1、掌握求数值导数和数值积分的方法

2、掌握代数方程数值求解的方法

3、掌握常微分方程数值求解的方法

二、实验内容

1、求函数在指定点的数值导数。

()

23

2

123,1,2,3 026

x x x

f x x x x

x

==

>>syms x

>>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f)

F=2*x^3

>>h=0.1

>>x=[0:h:4];

>>f=2*x^3;

>>[dy,dx]=diff_ctr(f,h,1);

>>y1=dy(dx==1)

y1=6.0000

>>y2=dy(dx==2)

y2=24.0000

>>y3=dy(dx==3) y3=54.0000

2、用数值方法求定积分。

（1）

210I π

=?的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)');

I=quadl(a,0,2*pi)

I =

6.7992 + 3.1526i

（2）()1

202ln 11x I dx x

+=+? b=inline('log(1+x)./(1+x.^2)');

I=quadl(b,0,1)

I =

0.2722

3、分别用3种不同的数值方法解线性方程组。

6525494133422139211

x y z u x y z u x y z u x y u +-+=-??-+-=??++-=??-+=?

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

b=[-4,13,1,11]';

x=A\b

x =

0.6667

-1.0000

1.5000

-0.0000

4、求非齐次线性方程组的通解。

123412341

234273635224

9472x x x x x x x x x x x x +++=??+++=??+++=?

A=[2 7 3 1;3 5 2 2 ;9 4 1 7]

A =

2 7

3 1

3 5 2 2

9 4 1 7

>> b=[6 4 2]'

b =

6

4

2

>> rank(A)

ans = 2

>> rank([A,b])

ans =

2

>> x=pinv(A)*b

x =

-0.0571

0.7429

0.3429 -0.1143

5、求代数方程的数值解。

（1）03sin 0 1.5x x x e x +-==在附近的根

function f=fz(x)

f=3*x-sin(x*pi/180)-exp(x)

fzero('fz',-1.5)

ans =

0.6287

（2）在给定的初值0001,1,1x y z ===下，求方程组的数值解。

2sin ln 703231050y x y z x z x y z ?++-=?+-+=??++-=?

function F=myfun (X) x=X(1);

y=X(2);

z=X(3);

F(1)=sin(x)+y^2+log(z)-7;

F(2)=3*x+2^y-3*z+1;

F(3)=x+y+z-5;

X=fsolve('myfun',[1,1,1],optimset('Display', 'off')) X = 0.2136 2.4359 2.3505

6、求函数在指定区间的极值。

（1）()3cos log x

x x x x f x e ++=在（0，1）内的最小值。 function f=f(x)

f=(x^3+cos(x)+x+x*log(x))/exp(x);

[x,fval]=fminbnd(f,0,1)

（2）()321,21121222410f x x x x x x x x =+-+在[0，

0]附近的最小值点和最小值。

function fxy(u)

x=u(1);

y=u(2);

f=2.*x.^2-6.*x.*y+y.^2;

[U,fmin]=fminsearch('fxy',[0,0])

7、求微分方程的数值解，并绘制解的曲线。

()()22'500000xd y dy y dx dx y y ?-+=???=??=???

function dy=vdp5(x,y)

dy=zeros(2,1)

dy(1)=y(2)

dy(2)=5*y(2)/x-y(1)/x

[X,Y]=ode23('vdp5',[0,100],[0,0]);

plot(X,Y(:,1),'-')

8、求微分方程组的数值解，并绘制解的曲线。()()()'1'213'3121230.5100,01,01y yy y y y y y y y y y ?=?=-??=-??===?

function dy=rigid(t,y)

dy=zeros(3,1);

dy(1)=y(2)*y(3);

dy(3)=-0.51*y(1)*y(2);

[T,Y]=ode45('rigid',[0,20],[0,1,1]);

plot(T,Y(:,1),'-',T,Y(:,2),'+',T,Y(:,3),'*')

用递推公式计算定积分(matlab版)

用递推公式计算定积分 实验目的： 1．充分理解不稳定的计算方法会造成误差的积累，在计算过程中会导致误差的迅速增加，从而使结果产生较大的误差。 2．在选择数值计算公式来进行近似计算时，应学会选用那些在计算过程中不会导致误差迅速增长的计算公式。 3．理解不稳定的计算公式造成误差积累的来源及具体过程； 4．掌握简单的matlab语言进行数值计算的方法。 实验题目： 对n=0,1,2,…,20，计算定积分： 实验原理： 由于y(n)= = – 在计算时有两种迭代方法，如下： 方法一： y(n)=– 5*y(n-1),n=1,2,3, (20) 取y(0)= = ln6-ln5 ≈ 0.182322 方法二： 利用递推公式：y(n-1)=-*y(n),n=20,19, (1) 而且，由 = * ≤≤* =

可取：y(20)≈*()≈0.008730. 实验容： 对算法一,程序代码如下： function [y,n]=funa() syms k n t; t=0.182322; n=0; y=zeros(1,20); y(1)=t; for k=2:20 y(k)=1/k-5*y(k-1); n=n+1; end y(1:6) y(7:11) 对算法二，程序代码如下： %计算定积分； %n--表示迭代次数； %y用来存储结果； function [y,n]=f(); syms k y_20;

y=zeros(21,1); n=1; y_20=(1/105+1/126)/2; y(21)=y_20; for k=21:-1:2 y(k-1)=1/(5*(k-1))-y(k)/5; n=n+1; end 实验结果： 由于计算过程中，前11个数字太小，后9个数字比较大，造成前面几个数字只显示0.0000的现象，所以先输出前6个，再输出7—11个，这样就能全部显示出来了。 算法一结果： [y,n]=funa %先显示一y(1)—y(6) ans = 0.1823 -0.4116 2.3914 -11.7069 58.7346

线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB实现综述廖淑芳20122090 数计学院12计算机科学与技术1班（职教本科）一、分析课题 随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。 直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。 迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。 二、研究课题-线性代数方程组数值解法 一、直接法 1、Gauss消元法 通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x向量。

1.1消元过程 1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。 11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ??? (1)(1)(1)(1)(1)11 121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()000 00 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ?? ? 步骤如下： 第一步：1 11 1,2,,i a i i n a -? +=第行第行 11121121222212 n n n n nn n a a a b a a a b a a a b ?? ? ? ? ??? 111211(2)(2)(2)2222 (2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ??? 第二步：(2)2 (2)222,3, ,i a i i n a -?+=第行第行 111211(2)(2)(2)2222 (2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ??? 11121311(2)(2)(2)(2)222322 (3)(3)(3)33 33(3)(3)(3)3 00000 n n n n nn n a a a a b a a a b a a b a a b ?? ? ? ? ? ? ?? ? 类似的做下去，我们有： 第k 步：() ()k ,1, ,k ik k kk a i i k n a -?+=+第行第行。 n －1步以后，我们可以得到变换后的矩阵为： 11121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()00000 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ?? ?

matlab数值微积分与方程数值求解

电子一班王申江 实验九数值微积分与方程数值求解 一、实验目的 1、掌握求数值导数和数值积分的方法 2、掌握代数方程数值求解的方法 3、掌握常微分方程数值求解的方法 二、实验内容 1、求函数在指定点的数值导数。 () 23 2 123,1,2,3 026 x x x f x x x x x == >>syms x >>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f) F=2*x^3 >>h=0.1 >>x=[0:h:4]; >>f=2*x^3; >>[dy,dx]=diff_ctr(f,h,1); >>y1=dy(dx==1) y1=6.0000 >>y2=dy(dx==2)

y2=24.0000 >>y3=dy(dx==3) y3=54.0000 2、用数值方法求定积分。 （1） 210I π =?的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)'); I=quadl(a,0,2*pi) I = 6.7992 + 3.1526i （2）()1 202ln 11x I dx x +=+? b=inline('log(1+x)./(1+x.^2)'); I=quadl(b,0,1) I = 0.2722 3、分别用3种不同的数值方法解线性方程组。 6525494133422139211 x y z u x y z u x y z u x y u +-+=-??-+-=??++-=??-+=? A=[6,5,-2,5;9,-1,4,-1;3,4,2,-2;3,-9,0,2]; b=[-4,13,1,11]'; x=A\b

matlab求定积分之实例说明

一、符号积分 符号积分由函数int来实现。该函数的一般调用格式为： int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分； int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分； int(s,v,a,b)：求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间[a,b]上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间[a,b]上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。 例： 求函数x^2+y^2+z^2的三重积分。内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下： >>syms x y z %定义符号变量 >>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式 F2 = 1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2 ^(3/4) %给出有理数解 >>VF2=vpa(F2) %给出默认精度的数值解 VF2 = 224.92153573331143159790710032805 二、数值积分 1.数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)?法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。 2.数值积分的实现方法 基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为： [I,n]=quad('fname',a,b,tol,trace) 基于变步长、牛顿－柯特斯(Newton-Cotes)法，MATLAB给出了quadl函数来求定积分。该函数的调用格式为： [I,n]=quadl('fname',a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。 例： 求函数'exp(-x*x)的定积分，积分下限为0，积分上限为1。 >>fun=inline('exp(-x.*x)','x'); %用内联函数定义被积函数fname

数值积分算法与MATLAB实现陈悦5133201讲解

东北大学秦皇岛分校 数值计算课程设计报告 数值积分算法及MATLAB实现 学院数学与统计学院 专业信息与计算科学 学号5133201 姓名陈悦 指导教师姜玉山张建波 成绩 教师评语： 指导教师签字： 2015年07月14日

1 绪论 数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算. 现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域. 1.1 数值积分介绍 数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系. 求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础. 构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(Approximations Using Parabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现

MATLAB数值微积分

4.1数值微积分 4.1.1近似数值极限及导数 Matlab 数值计算中，没有求极限指令，也没有求导指令，而是利用差分指令： 用一个简单矩阵表现diff和gradient指令计算方式。 差分： Dx=diff(X) 对向量： Dx=X(2:n)-X(1:n-1) 对矩阵： DX=X(2:n,:)-X(1:n-1,:) 长度小1. DIFF(X), for a vector X, is [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. DIFF(X), for a matrix X, is the matrix of row differences, （结果缺少一行） [X(2:n,:) - X(1:n-1,:)]. DIFF(X,N,DIM) is the Nth difference function along dimension DIM. If N >= size(X,DIM), DIFF returns an empty array （N阶差分）

梯度： FX=gradient(F) Fx(1)=Fx(2)-Fx(1); F=[1,2,3;4,5,6;7,8,9] Dx=diff(F) （按行） Dx_2=diff(F,1,2) （按列） [FX,FY]=gradient(F) Fx(1)=Fx(2)-Fx(1), Fx(end)=F(end)-F(end-1) FX与F维数相同。 [FX_2,FY_2]=gradient(F,0.5) %采样间隔0.5 即： Fx(1)=(Fx(2)-Fx(1))/2 F = 1 2 3 4 5 6 7 8 9 Dx = 3 3 3 3 3 3 Dx_2 = 1 1 1 1 1 1 FX = 1 1 1

实验一B Matlab基本操作与微积分计算

实验一Matlab基本操作与微积分计算 实验目的 1.进一步理解导数概念及其几何意义. 2.学习matlab的求导命令与求导法. 3.通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. 4.学习并掌握用matlab求不定积分、定积分、二重积分、曲线积分的方法. 5.学习matlab命令sum、symsum与int. 实验内容 一、变量 1、变量 MA TLAB中变量的命名规则是: （1）变量名必须是不含空格的单个词; （2）变量名区分大小写; （3）变量名最多不超过19个字符; （4）变量名必须以字母打头,之后可以是任意字母、数字或下划线,变量名中不允许使用标点符号. 1、创建简单的数组 x=[a b c d e f ]创建包含指定元素的行向量 x=first:step: last创建从first起,逐步加step计数,last结束的行向量, step缺省默认值为1 x=linspace(first,last,n）创建从first开始,到last结束,有n个元素的行向量 x=logspace(first,last,n）创建从first开始,到last结束,有n个元素的对数分隔行向量. 注:以空格或逗号分隔的元素指定的是不同列的元素,而以分号分隔的元素指定了不同行的元素. 2、数组元素的访问 （1）访问一个元素: x(i)表示访问数组x的第i个元素. （2）访问一块元素: x(a :b :c)表示访问数组x的从第a个元素开始,以步长为b到第c个元素（但

不超过c）,b可以为负数,b缺损时为1. （3）直接使用元素编址序号: x ([a b c d]) 表示提取数组x的第a、b、c、d个元素构成一个新的数组[x (a) x (b) x(c) x(d)]. 3、数组的运算 （1）标量-数组运算 数组对标量的加、减、乘、除、乘方是数组的每个元素对该标量施加相应的加、减、乘、除、乘方运算. 设:a=[a1,a2,…,an], c=标量, 则: a+c=[a1+c,a2+c,…,an+c] a .*c=[a1*c,a2*c,…,an*c] a ./c= [a1/c,a2/c,…,an/c](右除） a .\c= [c/a1,c/a2,…,c/an] (左除） a .^c= [a1^c,a2^c,…,an^c] c .^a= [c^a1,c^a2,…,c^an] （2）数组-数组运算 当两个数组有相同维数时,加、减、乘、除、幂运算可按元素对元素方式进行的,不同大小或维数的数组是不能进行运算的. 设:a=[a1,a2,…,an], b=[b1,b2,…,bn], 则: a +b= [a1+b1,a2+b2,…,an+bn] a .*b= [a1*b1,a2*b2,…,an*bn] a ./b= [a1/b1,a2/b2,…,an/bn] a .\b=[b1/a1,b2/a2,…,bn/an] a .^b=[a1^b1,a2^b2,…,an^bn] 三、矩阵 1、矩阵的建立 矩阵直接输入:从“[ ” 开始,元素之间用逗号“,”(或空格),行之间用分号“;”(或回车),用“ ]”结束. 特殊矩阵的建立: a=[ ] 产生一个空矩阵,当对一项操作无结果时,返回空矩阵,空矩阵的大小为零. b=zeros (m,n) 产生一个m行、n列的零矩阵 c=ones (m,n) 产生一个m行、n列的元素全为1的矩阵 d=eye (m,n) 产生一个m行、n列的单位矩阵 eye (n) %生成n维的单位向量 eye (size (A)) %生成与A同维的单位阵 2、矩阵中元素的操作 （1）矩阵A的第r行A（r,:） （2）矩阵A的第r列A（:,r） （3）依次提取矩阵A的每一列,将A拉伸为一个列向量A（:） （4）取矩阵A的第i1~i2行、第j1~j2列构成新矩阵:A(i1:i2, j1:j2) （5）以逆序提取矩阵A的第i1~i2行,构成新矩阵:A(i2:-1:i1,:） （6）以逆序提取矩阵A的第j1~j2列,构成新矩阵:A(:, j2:-1:j1 ） （7）删除A的第i1~i2行,构成新矩阵:A(i1:i2,:)=[ ] （8）删除A的第j1~j2列,构成新矩阵:A(:, j1:j2)=[ ] （9）将矩阵A和B拼接成新矩阵:[A B];[A;B] 3、矩阵的运算 （1）标量-矩阵运算同标量-数组运算. （2）矩阵-矩阵运算 a. 元素对元素的运算,同数组-数组运算.(A/B %A右除B; B\A%A左除B) b. 矩阵运算: 矩阵加法:A+B 矩阵乘法:A*B 方阵的行列式:det（A） 方阵的逆:inv（A）

用MATLAB算多元函数积分

用MATLAB 计算多元函数的积分 三重积分的计算最终是化成累次积分来完成的，因此只要能正确的得出各累次积分的积分限，便可在MA TLAB 中通过多次使用int 命令来求得计算结果。但三重积分的积分域Ω是一个三维空间区域，当其形状较复杂时，要确定各累次积分的积分限会遇到一定困难，此时，可以借助MATLAB 的三维绘图命令，先在屏幕上绘出Ω的三维立体图，然后执行命令 rotate3d on ↙ 便可拖动鼠标使Ω的图形在屏幕上作任意的三维旋转，并且可用下述命令将Ω的图形向三个坐标平面进行投影： view(0,0),向XOZ 平面投影； view(90,0),向YOZ 平面投影； view(0,90),向XOY 平面投影. 综合运用上述方法，一般应能正确得出各累次积分的积分限。 例11.6.1计算zdv Ω ???，其中Ω是由圆锥曲面222z x y =+与平面z=1围成的闭区域 解 首先用MA TLAB 来绘制Ω的三维图形，画圆锥曲面的命令可以是： syms x y z ↙ z=sqrt(x^2+y^2); ↙ ezsurf(z,[-1.5,1.5]) ↙ 画第二个曲面之前，为保持先画的图形不会被清除，需要执行命令 hold on ↙ 然后用下述命令就可以将平面z=1与圆锥面的图形画在一个图形窗口内： [x1,y1]=meshgrid(-1.5:1/4:1.5); ↙ z1=ones(size(x1)); ↙ surf(x1,y1,z1) ↙ 于是得到Ω的三维图形如图:

由该图很容易将原三重积分化成累次积分： 111zdv dy -Ω=???? 于是可用下述命令求解此三重积分： clear all ↙ syms x y z ↙ f=z; ↙ f1=int(f,z.,sqrt(x^2+ y^2),1); ↙ f2=int(f1,x,-sqrt(1- y^2), sqrt(1- y^2)); ↙ int(f2,y,-1,1) ↙ ans= 1/4*pi 计算结果为4 π 对于第一类曲线积分和第一类曲面积分，其计算都归结为求解特定形式的定积分和二重积分，因此可完全类似的使用int 命令进行计算，并可用diff 命令求解中间所需的各偏导数。 例11.6.2用MATLAB 求解教材例11.3.1 解 求解过程如下 syms a b t ↙ x=a*cos(t); ↙ y=a*sin(t); ↙ z=b*t; ↙ f=x^2 +y^2+z^2; ↙ xt=diff(x,t); ↙ yt=diff(y,t); ↙ zt=diff(z,t); ↙ int(f*sqrt(xt^2 +yt^2+zt^2),t,0,2*pi) ↙ ans= 2/3*( a^2 +b^2)^1/2*a^2*pi+8/3*( a^2 +b^2)^1/2*b^2*pi^3 对此结果可用factor 命令进行合并化简： factor （ans ） ans= 2/3*( a^2 +b^2)^1/2*pi*(3* a^2 +4*b^2*pi^2) 例11.6.3用MATLAB 求解教材例11.4.1 解 求解过程如下 syms x y z1 z2↙ f= x^2 +y^2; ↙ z1=sqrt(x^2 +y^2); ↙ z2=1; ↙ z1x=diff(z1,x); ↙ z1y=diff(z1,y); ↙ z2x=diff(z2,x); ↙ z2y=diff(z2,y); ↙

Matlab数值积分与数值微分

M a t l a b数值积分与数值微分 Matlab数值积分 1.一重数值积分的实现方法 变步长辛普森法、高斯-克朗罗德法、梯形积分法 1.1变步长辛普森法 Matlab提供了quad函数和quadl函数用于实现变步长 辛普森法求数值积分.调用格式为: [I,n]=Quad(@fname,a,b,tol,trace) [I,n]=Quadl(@fname,a,b,tol,trace) Fname是函数文件名，a,b分别为积分下限、积分上限； tol为精度控制，默认为1.0×10-6，trace控制是否展 开积分过程，若为0则不展开，非0则展开，默认不展开. 返回值I为积分数值；n为调用函数的次数. --------------------------------------------------------------------- 例如：求 ∫e0.5x sin(x+π )dx 3π 的值. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数

[I,n]=quad(@fesin,0,3*pi,1e-10) I= 0.9008 n= 365 --------------------------------------------------------------------- 例如：分别用quad函数和quadl函数求积分 ∫e0.5x sin(x+π 6 )dx 3π 的近似值，比较函数调用的次数. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6)); formatlong [I,n]=quadl(@fesin,0,3*pi,1e-10) I= n= 198 [I,n]=quad(@fesin,0,3*pi,1e-10) I= n= 365 --------------------------------------------------------------------- 可以发现quadl函数调用原函数的次数比quad少，并 且比quad函数求得的数值解更精确. 1.2高斯-克朗罗德法

matlab微积分基本运算

matlab 微积分基本运算 §1 解方程和方程组解 1. 线性方程组求解 对于方程 AX = B ,其中 A 是( m ×n )的矩阵有三种情形： 1)当n=m 且A 非奇异时,此方程为“恰定”方程组。 2)当 n > m 时,此方程为“超定”方程组。 3)当n

0.3188 两种方法所求方程组的解相同。 (2)MATLAB 解超定方程AX=B 的方法 对于方程 AX = B ,其中 A 是( m ×n )的矩阵， n > m ，如果A 列满秩，则此方程是没有精确解的。然而在实际工程应用中，求得其最小二乘解也是有意义的。基本解法有： 1)采用求伪逆运算解方程 x=pinv(A)*B 说明：此解为最小二乘解x=inv(A ’*A)*A*B,这里pinv(A) =inv(A ’*A)*A. 2)采用左除运算解方程 x=A\B 例2 “求伪逆”法和“左除”法求下列方程组的解 ??? ??=+=+=+1 221421 221 2121x x x x x x 命令如下： >> a=[1 2;2 4;2 2]; >> b=[1,1,1]'; >> xc=a\b %用左除运算解方程 运行得结果: xc = 0.4000 0.1000 >> xd=pinv(a)*b %用求伪逆运算解方程 运行得结果: xd = 0.4000 0.1000 >> a*xc-b %xc 是否满足方程ax=b 运行得结果: ans = -0.4000 0.2000 0.0000 可见xc 并不是方程的精确解。 (3) MATLAB 解欠定方程AX=B 的方法 欠定方程从理论上说是有无穷多个解的，如果利用求“伪逆”法和“左除”法来求解，只能得到其中一个解。基本方法： 1)采用求伪逆运算解方程 x=pinv(A)*B 2)采用左除运算解方程

数值积分的matlab实现

实验10 数值积分 实验目的： 1．了解数值积分的基本原理； 2．熟练掌握数值积分的MATLAB 实现； 3．会用数值积分方法解决一些实际问题。 实验内容： 积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分?1 0 d sin x x x 。这时我们一般考虑用数值方法计算其 近似值，称为数值积分。 10.1 数值微分简介 设函数()y f x =在* x 可导，则其导数为 h x f h x f x f h ) ()(lim )(**0* -+='→ （10.1） 如果函数()y f x =以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值 h x f h x f x f ) ()()(*** -+≈' （10.2） 表 10-1 一般的，步长h 越小，所得结果越精确。（10.2）式右端项的分子称为函数()y f x =在 *x 的差分，分母称为自变量在*x 的差分，所以右端项又称为差商。数值微分即用差商近似 代替微商。常用的差商公式为： 000()() ()2f x h f x h f x h +--'≈ （10.3） h y y y x f 243)(2 100-+-≈ ' （10.4）

h y y y x f n n n n 234)(12+-≈ '-- （10.5） 其误差均为2 ()O h ，称为统称三点公式。 10.2 数值微分的MATLAB 实现 MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x) 其中x 是n 维数组，dx 为1n -维数组[]21321,, ,n x x x x x x ---,这样基于两点的数值导 数可通过指令diff(x)/h 实现。对于三点公式，读者可参考例1的M 函数文件diff3.m 。 例1 用三点公式计算()y f x =在=x 1.0,1.2,1.4处的导数值，()f x 的值由下表给 解：建立三点公式的M 函数文件diff3.m 如下： function f=diff3(x,y) n=length(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h); for j=2:n-1 f(j)=(y(j+1)-y(j-1))/(2*h); end f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h); 在MATLAB 指令窗中输入指令： x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y) 运行得各点的导数值为：-0.2470，-0.2170，-0.1890，-0.1650，-0.0014。所以()y f x =在=x 1.0,1.2,1.4处的导数值分别为-0.2470，-0.1890和-0.0014。 对于高阶导数，MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如下： step1：对给定数据点（x,y ），利用指令pp=spline(x,y)，获得三次样条函数数据pp ，供后面ppval 等指令使用。其中，pp 是一个分段多项式所对应的行向量，它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。 step2：对于上面所求的数据向量pp ，利用指令[breaks,coefs,m,n]=unmkpp(pp)进行处理，生成几个有序的分段多项式pp 。 step3：对各个分段多项式pp 的系数，利用函数ppval 生成其相应导数分段多项式的系数，再利用指令mkpp 生成相应的导数分段多项式 step4：将待求点xx 代入此导数多项式，即得样条导数值。 上述过程可建立M 函数文件ppd.m 实现如下： function dy=ppd(pp) [breaks,coefs,m]=unmkpp(pp);

数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226 [摘要] 数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。 [关键字] 偏微分方程数值解matlab 数学物理方程的可视化 一：研究意义 在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。 二：数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题，将求解区域 {}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时 间[0,1]区间作n 等分，并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ 为空间和 时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2：差分格式的建立 0u u t x ??-=??………………………………（1） 设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)： ( )()G u u dxdt udt udx t x Γ??-=--??? ? 其中G Γ=?。于是可将（1）式写成积分守恒形式： ()0udt udx Γ --=? (2) 我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示

第七章 基于MATLAB的科学计算—数值微积分1

科学计算—理论、方法 及其基于MATLAB 的实现与分析 数值微积分 §1 数值微分 对于给定的函数()x f y =，如果 １、()x f y =的函数关系式比较复杂时； ２、 ()x f y =未知，而仅仅知道()x f y =在1+n 个相异点k x ， n k ,,1,0Λ=处 的函数值k y ； 则希望能用相对简单的计算方法，求得()x f y =导数的（近似）值。 基于上述考虑，选择的方法之一是利用函数()x f y =的插值多项式的导数作为函数()x f y =导数的近似值，例如Lagrange 插值多项式，由于 ()()()()∑==≈n k k k n x f x l x L x f 0 （１） 因而有 ()()x L x f n '≈' （２）

这里需要说明一点的是，尽管()x f 和()x L n 的函数值可能相差不多，但是它们的导数有可能相差很大，如下面的例子 例1： 考虑函数()2 2511x x f += 在区间[-1，1]的插值问题，取 区间[-1，1]的11个点 k x k ?+-=2.01，10,,1,0Λ=k ，作函数 ()2 2511x x f += 的10次插值多项式： ()18552.163597.1234338.3819095.4949417.220246810+-+-+-=x x x x x x L n 函数()x f 和插值多项式()x L n 的导数分别为 ()() 2 225150x x dx x df +-= ()x x x x x dx x dLn 7.334.4936.22883.39594.22093579-+-+-= 对函数()x f 和插值多项式()x L n 及其导数分别比较，结果如图所示： Derivative_Runge

详解Matlab求积分的各种方法

详解Matlab求积分的各种方法 一、符号积分由函数int来实现。 该函数的一般调用格式为： int(s)： 没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分；int(s,v)： 以v为自变量，对被积函数或符号表达式s求不定积分；int(s,v,a,b)： 求定积分运算。 a,b分别表示定积分的下限和上限。 该函数求被积函数在区间[a,b]上的定积分。 a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。 当函数f关于变量x在闭区间[a,b]上可积时，函数返回一个定积分结果。 当a,b中有一个是inf时，函数返回一个广义积分。 当a,b中有一个符号表达式时，函数返回一个符号函数。 例： 求函数x^2+y^2+z^2的三重积分。 内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下： >>syms x y z %定义符号变 量>>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式F2 =57/-

/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4) %给出有理数 解>>VF2=vpa(F2) %给出默认精度的数值解VF2 = 224.9 232805二、数值积分 1.数值积分基本原理求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)?法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。 它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]， i=1,2,…,n，其中x1=a，xn+1=b。 这样求定积分问题就分解为求和问题。 2.数值积分的实现方法基于变步长辛普生法，MATLAB给出了quad函数来求定积分。 该函数的调用格式为： [I,n]=quad('fname',a,b,tol,trace)基于变步长、牛顿－柯特斯(Newton-Cotes)法，MATLAB给出了quadl函数来求定积分。 该函数的调用格式为： [I,n]=quadl('fname',a,b,tol,trace)其中fname是被积函数名。 a和b分别是定积分的下限和上限。 tol用来控制积分精度，缺省时取tol= 0.0 01。 trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace= 0。

一维抛物线型方程数值解法(1)(附图及matlab程序)

一维抛物线偏微分方程数值解法（1） 解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法） Ut-Uxx=0, 00) U(x,0)=e^x, 0<=x<=1, U(0,t)=e^t,U(1,t)=e^(1+t), 00) %不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值 %m,n为x,t方向的网格数，例如（2-0）/0.01=200; %e为误差，p为精确解 u=zeros(n+1,m+1); x=0+(0:m)*h1; t=0+(0:n)*h2; for(i=1:n+1) u(i,1)=exp(t(i)); u(i,m+1)=exp(1+t(i)); end for(i=1:m+1) u(1,i)=exp(x(i)); end for(i=1:n+1) for(j=1:m+1) f(i,j)=0; end end r=h2/(h1*h1); %此处r=a*h2/(h1*h1）；a=1 要求r<=1/2差分格式才稳定for(i=1:n) for(j=2:m) u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j); end end for(i=1:n+1) for(j=1:m+1) p(i,j)=exp(x(j)+t(i)); e(i,j)=abs(u(i,j)-p(i,j)); end end

一维抛物线偏微分方程数值 解法(3)(附图及matlab程序)

一维抛物线偏微分方程数值解法（3） 上一篇参看一维抛物线偏微分方程数值解法（2）（附图及matlab程序） 解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法） Ut-Uxx=0, 00) U(x,0)=e^x, 0<=x<=1, U(0,t)=e^t,U(1,t)=e^(1+t), 0

r=h2/(h1*h1); for(i=1:n) %外循环，先固定每一时间层，每一时间层上解一线性方程组% a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2* (u(i+1,1)+u(i,1))+h2*f(i,j)... +(1-r)*u(i,2)+r/2*u(i,3); for(k=2:m-2) a(k)=-r/2;b(k)=1+r;c(k)=- r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)... *u(i,k+1)+r/2*u(i,k+2); %输入部分系数矩阵，为0的矩阵元素不输入% end a(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2* (u(i,m+1)+u(i+1,m+1)... )+r/2*u(i,m-1)+(1-r)*u(i,m); for(k=1:m-2) %开始解线性方程组消元过程 a(k+1)=-a(k+1)/b(k); b(k+1)=b(k+1)+a(k+1)*c(k); d(k+1)=d(k+1)+a(k+1)*d(k); end u(i+1,m)=d(m-1)/b(m-1); %回代过程% for(k=m-2:-1:1) u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k); end end for(i=1:n+1) for(j=1:m+1) p(i,j)=exp(x(j)+t(i)); %p为精确解 e(i,j)=abs(u(i,j)-p(i,j));%e为误差 end end [u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp; >> xlabel('x');ylabel('t');zlabel('e'); >> title('误差曲面')