matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

8.2.3 理查森外推法求导及其MATLAB程序

8.2.3  理查森外推法求导及其MATLAB程序

(一)一般形式的理查森外推法及其MATLAB程序

利用理查森外推法计算的近似值和误差估计的MATLAB程序

function [Dy,dy,jdw,n]=diffext1(fun,x0,jdwc,max1)

h=1;j=1; n=1;jdW=1;xdW=1; x1=x0+h;x2=x0-h;

Dy(1,1)=(feval(fun,x1)- feval(fun,x2))/(2*h);

while((jdW>jdwc)&(j<max1))

j;x1=x0+2^(-j)*h;x2=x0-2^(-j)*h;

Dy(j+1,1)=(feval(fun,x1)-feval(fun,x2))/(2^(1-j)*h);

for k=1:j

k;Dy(j+1,k+1)= Dy(j+1,k)+( Dy(j+1,k)- Dy(j,k))/(4^k-1);

end

jdW=abs(Dy(j+1,j+1)-Dy(j+1,j)); j=j+1;

end

[n,n]=size(Dy);jdw=abs(Dy(n,n)-Dy(n,n-1));

dy= Dy(n,n);

例 8.2.5  设,取精度jdwc=0.000 000 1,用理查森外推法计算的近似值和误差估计,将近似值与精确值4.46550187104484比较.

解 取最大计算次数max1=100,输入程序

>> x0=0.79;jdwc=0.0000001,max1=100;

[Dy,dy,jdw,n]=diffext1(@fun,x0,jdwc,max1),

wu=4.46550187104484-dy

运行后屏幕显示的近似值dy,dy与精确值的绝对误差wu,导数近似值的迭代矩阵Dy,jdW=|Dy(n,n)-Dy(n,n-1)|的值,最佳近似值dy的坐标n如下

Dy =

  Columns 1 through 4

0.95036708207779                  0                  0                  0

0.87447447334140  0.84917693709594                  0                  0

1.04543344913993  1.10241977440611  1.11930263022679                  0

3.33291848491782  4.09541349684379  4.29494641167297  4.34535345582291

4.16327772677060  4.44006414072153  4.46304085031338  4.46570901600608

4.38872901334021  4.46387944219674  4.46546712896176  4.46550564132126

4.44623223629059  4.46539997727405  4.46550134627921  4.46550188941123

  Columns 5 through 7

                  0                    0                  0

                  0                    0                  0

                  0                    0                  0

                  0                    0                  0

   4.46618099859504                   0                  0

   4.46550484377347   4.46550418282057                  0

   4.46550187469786   4.46550187179553   4.46550187123118

dy =             jdw =                 n=      wu =

 4.46550187123118  5.643530087695581e-010  7  -1.863398324530863e-010

 

 

(二)精度为的中心差商公式和误差估计及其MATLAB程序

用精度为的中心差商公式计算的近似值和误差估计的MATLAB程序

function [x0,yx,wuc]=zxcs4(h,x0,fi,M)

xi=[x0-2*h,x0-h,x0,x0+h,x0+2*h];

x1= xi(1); x2= xi(2); x3= xi(3); x4= xi(4);

x5= xi(5);y1=fi(1); y2=fi(2); y3=fi(3);

y4=fi(4); y5=fi(5);

yx=(8*y4-8*y2-y5+y1)/(12*h);

wuc=abs(h.^4*M/30);

用精度为的中心差商公式计算的近似值和误差估计的MATLAB程序

function [x,y,yx,wuc]=zxcs5(fun,h,x0, M)

x=zeros(1,5);y= zeros(1,5);

for k=1:5

x(k)=x0+(k-3).*h;

y(k)= feval(fun, x(k));

end

x;y;

yx=(8*y(4)-8*y(2)-y(5)+y(1))/(12*h);

wuc=abs(h.^4*M/30);

 

例 8.2.6  设.

(1)分别根据(8.7) 式和(8.24)式计算的近似值,并估计误差,取小数点后14位计算,其中步长分别取 M=872,.

(2)将(1)中计算的的近似值分别与精确值比较.

解 (1)计算的一阶导数的近似值和误差估计.

方法1 输入程序

>>x0=0.79;h1=0.1,M=872;

[x,y,yx1,wuc1]=zxcs5(@fun,h1,x0,M)

h2=0.01, [x,y,yx2,wuc2]=zxcs5(@fun,h2,x0, M),

h3=0.001,[ x,y,yx3,wuc3]=zxcs5(@fun,h3,x0, M),

h4=0.0001, [x,y,yx4,wuc4]=zxcs5(@fun,h4,x0, M),

运行后屏幕显示x=(-2h, -h, ,+h, +2h), 在x处函数值向量y,步长分别取时,根据 (8.24)式计算函数处的导数的近似值yx1, yx2, yx3,yx4和误差估计wuc1,wuc2,wuc3,wuc4如下

x =

  Columns 1 through 3

 0.59000000000000   0.69000000000000   0.79000000000000

    Columns 4 through 5

 0.89000000000000   0.99000000000000

y =

  Columns 1 through 3

-0.39855804055354   0.22803197210174   0.82491732446828

  Columns 4 through Column 5

0.97151370486625    0.38160930304430

yx1 =                     wuc1 =

   4.30640543209856           0.00290666666667

yx2 =                       wuc2 =

   4.46548476000396           2.906666666666667e-007

yx3 =                     wuc3 =

   4.46550186933213           2.906666666666667e-011

yx4 =                     wuc4 =

   4.46550187103480           2.906666666666667e-015

方法2编写计算的一阶导数的近似值和误差估计的MATLAB程序,并输入此程序

>> x=0.79;h=[0.1,0.01,0.001,0.0001]; M=872;x1=x+h;x2=x-h;

x3=x+2*h;x4=x-2*h; y1=sin(5.*x1.^2-21);y2=sin(5.*x2.^2-21);

y3=sin(5.*x3.^2-21); y4=sin(5.*x4.^2-21); yc2=(y1-y2)./(2*h),

wuc2=abs(h.^2*M/6),yc4=(8*y1-8*y2-y3+y4)./ (12*h),

wuc4=abs(h.^4*M/30),syms x,f=sin(5.*x.^2-21); dy=diff(f,x)

运行后屏幕显示步长分别取,M= 872,分别根据(8.7) 式和(8.24)式计算的近似值yc2,yc4及其误差估计值wuc2,wuc4,导函数dy如下

yc2 =

3.71740866382257  4.45760286140635  4.46542283857626  4.46550108070765

wuc2 =

1.45333333333333  0.01453333333333  0.00014533333333  0.00000145333333

yc4 =

4.30640543209856  4.46548476000396  4.46550186933213  4.46550187103480

wuc4 =

0.00290666666667  0.00000029066667  0.00000000002907  0.00000000000000

dy =

10*cos(5*x^2-21)*x

(2)计算的值.输入程序

>> x=0.79; dy =10*cos(5*x^2-21)*x,

wu2=abs(yc2-dy), wu4=abs(yc4-dy)

运行后屏幕显示的近似值与精确值的绝对误差wuc2,wuc4和的精确值dy如下

dy =

4.46550187104484

wuc2=

0.74809320722227  0.00789900963848  0.00007903246857  0.00000079033719

wuc4=

0.15909643894627  0.00001711104088  0.00000000171271  0.00000000001003

 

(三)变步长的中心差商公式及其MATLAB程序

用变步长的中心差商公式计算的近似值和误差估计的MATLAB程序

function [n,H,Dy,W]=difflim(fun,x0,h0,wu,max1)

Dy=zeros(1,max1); W=zeros(1, max1); H=zeros(1, max1);

h=h0;H(1)=h; E(1)=0; x1=x0+h;x2=x0-h;

h1=h/10;x3=x0+h1;x4=x0-h1;

Dy(1)=(feval(fun,x1)- feval(fun,x2))./(2*h);

Dy(2)=(feval(fun,x3)- feval(fun,x4))./(2*h1);

W(1)=abs(Dy(2)- Dy(1));

k=1;

while( (W(k)>wu)&(k<max1) )

h=h/(10^(k)) ; H(k+1)=h; x1=x0+H(k+1); x2=x0-H(k+1);

Dy(k+1)=(feval(fun,x1)- feval(fun,x2))./(2* H(k+1));

W(k+1)=abs(Dy(k+1)- Dy(k));

k=k+1;

end

n=length(Dy(1:k))-2;

Dy= Dy(2:k); W=W(2:k); H=H(2:k);

 

例 8.2.7  设.取初始步长h0=0.2,精度wu=0.00001, 用变步长的中心差商公式计算的近似值和误差估计,与精确值4.46550187104484比较.

解  输入程序

>> x0=0.79;wu=0.00001;max1=100;h0=0.2;

[n,H,Dy,W]= difflim(@fun,x0,h0,wu,max1)

jdwc=4.46550187104484- Dy

运行后屏幕显示迭代次数n, 变步长数组H,用变步长的中心差商公式计算的近似值Dy,误差估计值W,jdwc,精确值与近似值Dy的差如下

n =

     2

H =

   0.02000000000000   0.00020000000000   0.00000020000000

Dy =

   4.43395716561354   4.46549870972618   4.46550187410688

W =

   2.48353880661895   0.03154154411264   0.00000316438070

jdwc =

   0.03154470543130   0.00000316131866  -0.00000000306204

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3080203.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1