第三章微积分问题的计算机求解

第三章微积分问题的计算机求解内容摘要：

sin(i*pi*x/L)。 end if a+b, F=subs(F,x,xLa)。 end ％换回变量区域 例 : syms x。 f=x*(xpi)*(x2*pi)。 [A,B,F]=fseries(f,x,6,0,2*pi) A = [ 0, 0, 0, 0, 0, 0, 0] B = [ 12, 3/2, 4/9, 3/16, 12/125, 1/18] F = 12*sin(x)+3/2*sin(2*x)+4/9*sin(3*x)+3/16*sin(4*x)+12/125*sin(5*x)+1/18*sin(6*x) 例 : syms x。 f=abs(x)/x。 % 定义方波信号 xx=[pi:pi/200:pi]。 xx=xx(xx~=0)。 xx=sort([xx,eps,eps])。 % 剔除零点 yy=subs(f,x,xx)。 plot(xx,yy,39。 r.39。 ), hold on % 绘制出理论值并保持坐标系 for n=2:20 [a,b,f1]=fseries(f,x,n), y1=subs(f1,x,xx)。 plot(xx,y1) end a = [ 0, 0, 0] b = [ 4/pi, 0] f1 = 4/pi*sin(x) a = [ 0, 0, 0, 0] b = [ 4/pi, 0, 4/3/pi] f1 = 4/pi*sin(x)+4/3/pi*sin(3*x) …… 级数求和的计算 • 是在符号工具箱中提供的 例 :计算 format long。 sum(2.^[0:63]) %数值计算 ans = +019 sum(sym(2).^[0:200]) % 或 syms k。 symsum(2^k,0,200) ％把 2定义为符号量可使计算更精确 ans = 3213876088517980551083924184682325205044405987565585670602751 syms k。 symsum(2^k,0,200) ans = 3213876088517980551083924184682325205044405987565585670602751 例 :试求解无穷级数的和 syms n。 s=symsum(1/((3*n2)*(3*n+1)),n,1,inf) %采用符号运算工具箱 s = 1/3 m=1:10000000。 s1=sum(1./((3*m2).*(3*m+1)))。 %数值计算方法 ,双精度有效位 16,“大数吃小数” ,无法精确 format long。 s1 % 以长型方式显示得出的结果 s1 = 例 :求解 syms n x s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf)。 simple(s1) % 对结果进行化简， MATLAB 及以前版本因本身 bug 化简很麻烦 ans = log((((2*x+1)^2)^(1/2)+1)/(((2*x+1)^2)^(1/2)1)) %实际应为 log((x+1)/x) 例 :求 syms m n。 limit(symsum(1/m,m,1,n)log(n),n,inf) ans = eulergamma vpa(ans, 70) % 显示 70 位有效数字 ans = .5772156649015328606065120900824024310421593359399235988057672348848677 数值微分 39。 0( ) ( ) ( ) l imhf x h f xfxh差商型求导公式 由导数定义39。 39。 39。 1 ( ) ( ) ( )2 ( ) ( ) ( )3 ( ) ( ) ( )2f x h f xfxhf x f x hfxhf x h f x hfxh  （）向前差商公式（ ）向后差商公式（）中心差商公式 （中点方法 ） xh x x+h B C A T f(x) 数值微分算法 • 向前差商公式： • 向后差商公式 两种中心公式： 39。 2 39。 39。 3 39。 39。 39。 439。 2 39。 39。 3 39。 39。 39。 4239。 39。 39。 39。 ( ) ( ) ( ) / 2 ! ( ) / 3! ( )()2( ) ( ) ( ) / 2 ! ( ) / 3! ( )2( ) ( )3!f x t f x t f x t f o tfxtf x t f x t f x t f o tttf x f               中心差分方法及其 MATLAB 实现 function [dy,dx]=diff_ctr(y, Dt, n) yx1=[y 0 0 0 0 0]。 yx2=[0 y 0 0 0 0]。 yx3=[0 0 y 0 0 0]。 yx4=[0 0 0 y 0 0]。 yx5=[0 0 0 0 y 0]。 yx6=[0 0 0 0 0 y]。 switch n case 1 dy = (diff(yx1)+7*diff(yx2)+7*diff(yx3) … diff(yx4))/(12*Dt)。 L0=3。 case 2 dy=(diff(yx1)+15*diff(yx2) 15*diff(yx3)… +diff(yx4))/(12*Dt^2)。 L0=3。 ％数值计算 diff(X)表示数组 X相邻两数的差 case 3 dy=(diff(yx1)+7*diff(yx2)6*diff(yx3)6*diff(yx4)+... 7*diff(yx5)diff(yx6))/(8*Dt^3)。 L0=5。 case 4 dy = (diff(yx1)+11*diff(yx2)28*diff(yx3)+28*… diff(yx4)11*diff(yx5)+diff(yx6))/(6*Dt^4)。 L0=5。 end dy=dy(L0+1:endL0)。 dx=([1:length(dy)]+L02(n2))*Dt。 调用格式： y为 等距实测数据， dy为得出的导数向量， dx为相应的自变量向量， dy、 dx的数据比 y短。 [ , ] _ ( , , )yxd d d iff c tr y t n• 例： 求导数的解析解 ,再 用数值微分求取原函数的 1～4 阶导数，并和解析解比较精度。 h=。 x=0:h:pi。 syms x1。 y=sin(x1)/(x1^2+4*x1+3)。 % 求各阶导数的解析解与对照数据 yy1=diff(y)。 f1=subs(yy1,x1,x)。 yy2=diff(yy1)。 f2=subs(yy2,x1,x)。 yy3=diff(yy2)。 f3=subs(yy3,x1,x)。 yy4=diff(yy3)。 f4=subs(yy4,x1,x)。 y=sin(x)./(x.^2+4*x+3)。 % 生成已知数据点 [y1,dx1]=diff_ctr(y,h,1)。 subplot(221),plot(x,f1,dx1,y1,39。 :39。 )。 [y2,dx2]=diff_ctr(y,h,2)。 subplot(222),plot(x,f2,dx2,y2,39。 :39。 ) [y3,dx3]=diff_ctr(y,h,3)。 subplot(223),plot(x,f3,dx3,y3,39。 :39。 )。 [y4,dx4]=diff_ctr(y,h,4)。 subplot(224),plot(x,f4,dx4,y4,39。 :39。 ) 求最大相对误差： norm((y4… f4(4:60))./f4(4:60)) ans = 用插值、拟合多项式的求导数 • 基本思想：当已知函数在一些离散点上的函数值时，该函数可用插值或拟合多项式来近似，然后对多项式进行微分求得导数。 • 选取 x=0附近的少量点 • 进行多项式拟合或插值 • g(x)在 x=0处的 k阶导数为 ( , ) , 1 , 2 , , 1iix y i n11 2 1()nnnng x c x c x c x c    () 1( 0 ) ! 0 , 1 , 2 , ,k nkg c k k n• 通过坐标变换用上述方法计算任意 x点处的导数值 • 令 • 将 g(x)写成 z的表达式 • 导数为 • 可直接用 拟合节点 得到系数 d=polyfit(xa,y,length(xd)1) z x a11 2 1( ) ( )nnnng x g z d z d z d z d     ( ) ( )1( ) ( 0 ) ! 0 , 1 , ,kknkg a g d k k n  ()gz ( , )iix a y id• 例：数据集合如下： xd: 0 yd: 计算 x=a=。 xd=[ 0 ]。 yd=[ ]。 a=。 L=length(xd)。 d=polyfit(xda,yd,L1)。 fact=[1]。 for k。