科学计算科学计算 (35).pdf

《科学计算科学计算 (35).pdf》由会员分享，可在线阅读，更多相关《科学计算科学计算 (35).pdf（21页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、Chapter 6 Numerical Calculus and Solutions Chapter 6 Numerical Calculus and Solutions of Equationsof Equations6.1 Numerical Differentiation and Integration6.2 Solve System of Linear Equations6.3 Application Examples of System of Linear Equations6.4 Solve Nonlinear Equation and Find Extreme Value of

2、Function6.5 Numerical Solution for Ordinary Differential Equations6.6 Application Examples of Ordinary Differential Equations6.16.1 Numerical Differentiation and Numerical Differentiation and IntegrationIntegrationNumerical differentiationNumerical integration=+0000()lim()()hfxf xhf xh=0000()lim()()

3、hfxf xf xhh=+0000()lim(/2)(/2)hfxf xhf xhh1 1Numerical DifferentiationNumerical Differentiation（1）Numerical difference and difference quotient In calculus,the derivative of any function f(x)at x0is defined by the limit:=+f xf xhf x000()()()=f xf xf xh000()()()=+f xf xhf xh000()(/2)(/2)If the limit a

4、s h approaches zero in the definition of limit is removed,the forward difference,the backward difference and the center difference with h as the step size at x0are obtained.The forward difference:The backward difference:The center difference:+000()()()fxf xhf xh000()()()fxf xf xhh+000()(/2)(/2)fxf x

5、hf xhhThe differentiation of the function f(x)at x0is close to the difference of function at this point,and the derivative of the f(x)at x0is close to the difference quotient of the function at this point.When the step size h is small enough,the forward difference quotient,the backward difference qu

6、otient and the central difference quotient with h as the step size at x0are obtained.The forward difference quotient:The backward difference quotient:The central difference quotient:（2）Implementation of numerical differentiationThe diff function in MATLAB can find forward difference,which has three

7、forms of syntax.dx=diff(x):Compute the first forward difference of vector x,dx.dx(i)=x(i+1)-x(i),i=1,2,n-1.dx=diff(x,n):Compute the n-order forward difference of vector x.For example,diff(x,2)=diff(diff(x).dx=diff(A,n,dim):Compute n-order difference of matrix A.Compute difference bycolumn when dim e

8、quals 1,and compute difference by row when dim equals 2.The diff function computes the difference between vector elements.Therefore,the number ofdifference vector elements is one fewer than that of original vector.Equally,if it is a matrix,thedifference matrix has one row or one column fewer than th

9、e original matrix.After computingthe difference,the difference quotient of f(x)at a certain point can be used as the approximatevalue of its derivative.Example 1:Set f(x)=sinx.Create random samplings in the range 0,2,and approximate derivatives of the function f(x),and compare it with the theoretica

10、l value of cosx.x=0,sort(2*pi*rand(1,5000),2*pi;y=sin(x);f1=diff(y)./diff(x);f2=cos(x(1:end-1);plot(x(1:end-1),f1,x(1:end-1),f2);d=norm(f1-f2)d=0.0433In some circumstances,it has some difficulties when using Newton-Leibniz formula.For example,when the original function of the integrand cannot be exp

11、ressed by elementary function or the integrand is presented in tabular form,it is necessary to use numerical methods to approximate definite integral.()d()()f xxF bF aab=2 2Numerical Numerical I Integrationntegration（1）Fundamental theorem of numerical integrationIn advanced mathematics,we could use

12、fundamental theorem of calculus to find definite integral.It can be computed with Newton-Leibniz formula as long as the original function of the integrand f(x)is found:There are numerical methods to find definite integral,such as trapezoid method,Simpson method,Gauss quadrature formula,etc.Their bas

13、ic idea is to divide the integral interval a,b into n subintervals xi,xi+1,in which i varied between 1 and n,and x1=a,xn+1=b.In this way,we turn a definite integral problem into a summation problem:In every subinterval,the value of definite integral can be approximated.So we can also avoid the diffi

14、culty of finding the original function with Newton-Leibniz formula.()d()d1 1 Sf xxf xxiixxinab=+=Based on the adaptive Simpson methodI,n=quad(filename,a,b,tol,trace)Based on the adaptive Gauss-Lobatto methodI,n=quadl(filename,a,b,tol,trace)In which the filename is the name of the integrand;a and b i

15、s the lower and upper limit respectively of definite integral.Integral limit a,b must be finite,and cant be infinite.Tol is used to control the precision of integral.Set tol equal to 10-6by default.Trace is used to determine if the integral process is presented.If it is nonzero,the integral process

16、is presented,If not,the process will not be presented.Set trace equal to 0 by default.Return argument I,namely the value of definite integral,and n represents the calling times of the integrand.（2）Implementation of numerical integration+2 0 141xdx format long f=(x)4./(1+x.2);I,n=quad(f,0,1,1e-8)I=3.

17、141592653733437n=61 I,n=quadl(f,0,1,1e-8)I=3.141592653589806n=48(atan(1)-atan(0)*4ans=3.141592653589793 format shortExample 2:Use the quad function and the quadlfunction to approximate definite integral.And compare calling times of the integrand with the same integral precision.Based on global adapt

18、ive quadrature methodI=integral(filename,a,b)In which,I is the integral we got,nd filename is the integrand;a and b is the lower and upper limit respectively of the definite integral.Integral limit can be infinite.11ln21=eIxxdxExample 3:Find the definite integral.File of the integrand,fe.m：function

19、f=fe(x)f=1./(x.*sqrt(1-log(x).2);I=integral(fe,1,exp(1)I=1.5708The quadgk function which is based on Gauss-Kronrod methodI,err=quadgk(filename,a,b)The err argument will return approximate error range.Other arguments have the same meaning and usage as the quad functions.The upper and lower limits of

20、the integral can be infinite or can be a complex.If they are complex,then the quadgk function will find integral in the complex plane.File of integrand fsx.m:function f=fsx(x)f=sin(1./x)./x.2;xx1sin1xd2 2+I=quadgk(fsx,2/pi,+Inf)I=1.0000Example 4:Find the definite integral.Trapezoidal integration fun

21、ction trapzgiven(xi,yi)(i=1,2,n),and a=x1x2 x=1:6;y=6,8,11,7,5,2;plot(x,y,-ko);grid on axis(1,6,0,11);I1=trapz(x,y)I1=35 I2=sum(diff(x).*(y(1:end-1)+y(2:end)/2)I2=35（3）Numerical solution of multiple definite integralsI=integral3(filename,a,b,c,d,e,f)I=triplequad(filename,a,b,c,d,e,f,tol)f x yx ycadb

22、(,)d d f x y zx y zecafdb(,)d d dNumerically evaluate double integral:I=integral2(filename,a,b,c,d)I=quad2d(filename,a,b,c,d)I=dblquad(filename,a,b,c,d,tol)Numerically evaluate triple integral:f1=(x,y)exp(-x.2/2).*sin(x.2+y);I1=quad2d(f1,-2,2,-1,1)I1=1.574 f2=(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x);I2=integral3(f2,0,pi,0,pi,0,1)I2=1.7328xxyxy2esin()dd/22 2 2 1 1+Example 6:Find double integral and triple integral respectively.xzx y zz yx4ed d d000122