中科院的matlab课件,中科院的matlab课件
PPT內(nèi)容
這是中科院的matlab課件，關(guān)于微積分問(wèn)題的計(jì)算機(jī)求解，包括了微積分問(wèn)題的解析解，函數(shù)的級(jí)數(shù)展開與級(jí)數(shù)求和問(wèn)題求解，數(shù)值微分，數(shù)值積分問(wèn)題，曲線積分與曲面積分的計(jì)算等內(nèi)容，歡迎點(diǎn)擊下載。
第三章 微積分問(wèn)題的計(jì)算機(jī)求解
微積分問(wèn)題的解析解
函數(shù)的級(jí)數(shù)展開與級(jí)數(shù)求和問(wèn)題求解
數(shù)值微分
數(shù)值積分問(wèn)題
曲線積分與曲面積分的計(jì)算
3.1 微積分問(wèn)題的解析解 3.1.1 極限問(wèn)題的解析解
單變量函數(shù)的極限
格式1： L= limit( fun, x, x0)
格式2： L= limit( fun, x, x0, ‘left’ 或 ‘right’)
例: 試求解極限問(wèn)題
>> syms x a b;
>> f=x*(1+a/x)^x*sin(b/x);
>> L=limit(f,x,inf)
L =
exp(a)*b
例:求解單邊極限問(wèn)題
>> syms x;
>> limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right')
ans =
12
在(-0.1,0.1)區(qū)間繪制出函數(shù)曲線：
>> x=-0.1:0.001:0.1;
>> y=(exp(x.^3)-1)./(1-cos(sqrt(x-sin(x))));
Warning: Divide by zero.
(Type "warning off
MATLAB:
divideByZero" to
suppress this warning.)??
>> plot(x,y,'-',[0],
[12],'o')
多變量函數(shù)的極限：
格式： L1=limit(limit(f,x,x0),y,y0)
或? L1=limit(limit(f,y,y0), x,x0)
如果x0 或y0不是確定的值，而是另一個(gè)變量的函數(shù)，如x->g(y)，則上述的極限求取順序不能交換。
例：求出二元函數(shù)極限值
>> syms x y a;
>> f=exp(-1/(y^2+x^2)) … *sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);
>> L=limit(limit(f,x,1/sqrt(y)),y,inf)
L =
exp(a^2)
3.1.2 函數(shù)導(dǎo)數(shù)的解析解
函數(shù)的導(dǎo)數(shù)和高階導(dǎo)數(shù)
格式： y=diff(fun,x)?????? %求導(dǎo)數(shù)
y= diff(fun,x,n)?????? %求n階導(dǎo)數(shù)
例：
一階導(dǎo)數(shù)：
>> syms x; f=sin(x)/(x^2+4*x+3);
>> f1=diff(f); pretty(f1)
cos(x)??????????????? sin(x) (2 x + 4)
---------------? -? -------------------
2???????????????????????? 2?????????????? 2
x? + 4 x + 3?????? (x? + 4 x + 3)
原函數(shù)及一階導(dǎo)數(shù)圖：
>> x1=0:.01:5;
>> y=subs(f, x, x1);
>> y1=subs(f1, x, x1);
>> plot(x1,y,x1,y1,‘:’)
更高階導(dǎo)數(shù)：
>> tic, diff(f,x,100); toc
elapsed_time =
4.6860
原函數(shù)4階導(dǎo)數(shù)
>> f4=diff(f,x,4);? pretty(f4)
2
sin(x)??????? cos(x) (2 x + 4)??????? sin(x) (2 x + 4)
------------ + 4 ------------------- - 12 -----------------
2????????????????????? 2??????????????? 2???????????? 2??????????????? 3
x? + 4 x + 3???? (x? + 4 x + 3)??????????? (x? + 4 x + 3)
3
sin(x)???????????? cos(x) (2 x + 4)?????? cos(x) (2 x + 4)
+ 12 --------------- - 24 ----------------- + 48 ----------------
2???????????????? 2????????? 2?????????????? 4????????????? 2?????????????? 3
(x? + 4 x + 3)???????? (x? + 4 x + 3)??????????? (x? + 4 x + 3)
4???????????????????????????? 2
sin(x) (2 x + 4)?????? sin(x) (2 x + 4)?????????? sin(x)
+ 24 ----------------- - 72 ----------------- + 24 ---------------
2??????????????? 5???????????? 2?????????????? 4?????????? 2?????????????? 3
(x? + 4 x + 3)??????????? (x? + 4 x + 3)???????? (x? + 4 x + 3)
多元函數(shù)的偏導(dǎo)：
格式：? f=diff(diff(f,x,m),y,n)
或?? f=diff(diff(f,y,n),x,m)
例：?????????????????????????????? 求其偏導(dǎo)數(shù)并用圖表示。
>> syms x y? z=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> zx=simple(diff(z,x))
zx =
-exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)
>> zy=diff(z,y)
zy =
(x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-x*y)
直接繪制三維曲面
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,z), axis([-3 3 -2 2 -0.7 1.5])
>> contour(x,y,z,30), hold on?? % 繪制等值線
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);??? % 偏導(dǎo)的數(shù)值解
>> quiver(x,y,zx,zy)? % 繪制引力線
例
>> syms x y z; f=sin(x^2*y)*exp(-x^2*y-z^2);
>> df=diff(diff(diff(f,x,2),y),z); df=simple(df);
>> pretty(df)
2?????? 2?????????? 2??????????????????? 2????????? 2
-4 z exp(-x? y - z ) (cos(x? y) - 10 cos(x? y) y x? + 4
2????? 4?? 2??????????? 2????? 4?? 2???????? 2
sin(x? y) x? y+ 4 cos(x? y) x? y? - sin(x? y))
多元函數(shù)的Jacobi矩陣:
格式：J=jacobian(Y,X)
其中，X是自變量構(gòu)成的向量，Y是由各個(gè)函數(shù)構(gòu)成的向量。
例：
試推導(dǎo)其 Jacobi 矩陣
>> syms r theta phi;
>> x=r*sin(theta)*cos(phi);
>> y=r*sin(theta)*sin(phi);
>> z=r*cos(theta);
>> J=jacobian([x; y; z],[r theta phi])
J =
[??? sin(theta)*cos(phi),? r*cos(theta)*cos(phi), -r*sin(theta)*sin(phi)]
[??? sin(theta)*sin(phi),? r*cos(theta)*sin(phi),? r*sin(theta)*cos(phi)]
[???????????? cos(theta),????????? -r*sin(theta),??????????????????????? 0??????????????????? ]
隱函數(shù)的偏導(dǎo)數(shù)：
格式：F=-diff(f,xj)/diff(f,xi)
例：
>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> pretty(-simple(diff(f,x)/diff(f,y)))
3????? 2????????? 2
-2 x + 2 + 2 x? + x? y - 4 x? - 2 x y
-? -----------------------------------------
x (x - 2) (2 y + x)
3.1.3 積分問(wèn)題的解析解
不定積分的推導(dǎo)：
格式： F=int(fun,x)
例：
用diff() 函數(shù)求其一階導(dǎo)數(shù)，再積分，檢驗(yàn)是否可以得出一致的結(jié)果。
>> syms x; y=sin(x)/(x^2+4*x+3); y1=diff(y);
>> y0=int(y1); pretty(y0)??? % 對(duì)導(dǎo)數(shù)積分
sin(x)?????? sin(x)
- 1/2 ------ + 1/2 ------
x + 3??????? x + 1
對(duì)原函數(shù)求4 階導(dǎo)數(shù)，再對(duì)結(jié)果進(jìn)行4次積分
>> y4=diff(y,4);
>> y0=int(int(int(int(y4))));
>> pretty(simple(y0))
sin(x)
------------
2
x? + 4 x + 3
例:說(shuō)明
>> syms a x; f=simple(int(x^3*cos(a*x)^2,x))
f = 1/16*(4*a^3*x^3*sin(2*a*x)+2*a^4 *x^4+6*a^2*x^2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a^4
>> f1=x^4/8+(x^3/(4*a)-3*x/(8*a^3))*sin(2*a*x)+...
(3*x^2/(8*a^2)-3/(16*a^4))*cos(2*a*x);
>> simple(f-f1)? % 求兩個(gè)結(jié)果的差
ans =
-3/16/a^4
定積分與無(wú)窮積分計(jì)算:
格式：? I=int(f,x,a,b)
格式： I=int(f,x,a,inf)
例：
>> syms x; I1=int(exp(-x^2/2),x,0,1.5)? ％無(wú)解
I1 =
1/2*erf(3/4*2^(1/2))*2^(1/2)*pi^(1/2)
>> vpa(I1,70)
ans =
1.085853317666016569702419076542265042534236293532156326729917229308528
>> I2=int(exp(-x^2/2),x,0,inf)
I2 =
1/2*2^(1/2)*pi^(1/2)
多重積分問(wèn)題的MATLAB求解
例：
>> syms x y z; f0=-4*z*exp(-x^2*y-z^2)*(cos(x^2*y)-10*cos(x^2*y)*y*x^2+...
4*sin(x^2*y)*x^4*y^2+4*cos(x^2*y)*x^4*y^2-sin(x^2*y));
>> f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);
>> f1=simple(int(f1,x))
f1 =
exp(-x^2*y-z^2)*sin(x^2*y)
>> f2=int(f0,z); f2=int(f2,x); f2=int(f2,x);
>> f2=simple(int(f2,y))
f2 =2*exp(-x^2*y-z^2)*tan(1/2*x^2*y)/(1+tan(1/2*x^2*y)^2) ???
f2=sin(x^2*y)/exp(y*x^2 + z^2)
>> simple(f1-f2)
ans =
0
順序的改變使化簡(jiǎn)結(jié)果不同于原函數(shù)，但其誤差為0，表明二者實(shí)際完全一致。這是由于積分順序不同，得不出實(shí)際的最簡(jiǎn)形式。
例：
>> syms x y z
>> int(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi)
ans =
(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2))*pi^2*hypergeom([1],[2],-pi^2) ???
Ei(n,z)為指數(shù)積分，無(wú)解析解，但可求其數(shù)值解：
>> vpa(ans,60)
ans =
3.10807940208541272283461464767138521019142306317021863483588
3.2 函數(shù)的級(jí)數(shù)展開與??? 級(jí)數(shù)求和問(wèn)題求解
3.2.1 Taylor 冪級(jí)數(shù)展開
3.2.2 Fourier 級(jí)數(shù)展開
3.2.3 級(jí)數(shù)求和的計(jì)算
3.2.1 Taylor 冪級(jí)數(shù)展開? 3.2.1.1 單變量函數(shù)的 Taylor? 冪級(jí)數(shù)展開
例:
>> syms x; f=sin(x)/(x^2+4*x+3);
>> y1=taylor(f,x,9); pretty(y1)
23?????? 34???????????? 4087?????? 3067???????? 515273????????? 386459
1/3 x - 4/9 x2? + -- x3?? - ---- x4? +? ------x5? - ------ x6? +---------- x7? - --------- x8
54???????? 81??????????? 9720????? 7290????????? 1224720??????? 918540
>> taylor(f,x,9,2)
ans =
1/15*sin(2)+(1/15*cos(2)-8/225*sin(2))*(x-2)+ (-127/6750*sin(2)-8/225*cos(2))*(x-2)^2 +(23/6750*cos(2)+628/50625*sin(2))*(x-2)^3 +(-15697/6075000*sin(2)+28/50625*cos(2))*(x-2)^4 +(203/6075000*cos(2)+6277/11390625*sin(2))*(x-2)^5 +(-585671/2733750000*sin(2)-623/11390625*cos(2))*(x-2)^6 +(262453/19136250000*cos(2)+397361/5125781250*sin(2))*(x-2)^7 +(-875225059/34445250000000*sin(2)-131623/35880468750*cos(2))*(x-2)^8
>> syms a; taylor(f,x,5,a)? % 結(jié)果較冗長(zhǎng)，顯示從略
ans =
sin(a)/(a^2+3+4*a) +(cos(a)-sin(a)/(a^2+3+4*a)*(4+2*a))/(a^2+3+4*a)*(x-a) +(-sin(a)/(a^2+3+4*a)-1/2*sin(a)-(cos(a)*a^2+3*cos(a)+4*cos(a)*a-4*sin(a)-2*sin(a)*a)/(a^2+3+4*a)^2*(4+2*a))/(a^2+3+4*a)*(x-a)^2+…
例:對(duì)y=sinx進(jìn)行Taylor冪級(jí)數(shù)展開,并觀察不同階次的近似效果。
>> x0=-2*pi:0.01:2*pi; y0=sin(x0); syms x; y=sin(x);
>>? plot(x0,y0,'r-.'), axis([-2*pi,2*pi,-1.5,1.5]); hold on
>> for n=[8:2:16]
p=taylor(y,x,n), y1=subs(p,x,x0); line(x0,y1)?? end
p =
x-1/6*x^3+1/120*x^5-1/5040*x^7
p =
x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9
p =
x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11
p =
x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13
p =
x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13-1/1307674368000*x^15
3.2.1.2 多變量函數(shù)的Taylor??????? 冪級(jí)數(shù)展開
多變量函數(shù)????????????????????????? 在
的Taylor冪級(jí)數(shù)的展開
例：？？？
>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> F=maple('mtaylor',f,'[x,y]',8)
F =
mtaylor((x^2-2*x)*exp(-x^2-y^2-x*y),[x, y],8)
>> maple(‘readlib(mtaylor)’);％讀庫(kù)，把函數(shù)調(diào)入內(nèi)存
>> F=maple('mtaylor',f,'[x,y]',8)
F =
-2*x+x^2+2*x^3-x^4-x^5+1/2*x^6+1/3*x^7+2*y*x^2+2*y^2*x-y*x^3-y^2*x^2-2*y*x^4-3*y^2*x^3-2*y^3*x^2-y^4*x+y*x^5+3/2*y^2*x^4+y^3*x^3+1/2*y^4*x^2+y*x^6+2*y^2*x^5+7/3*y^3*x^4+2*y^4*x^3+y^5*x^2+1/3*y^6*x
>> syms a; F=maple('mtaylor',f,'[x=1,y=a]',3);
>> F=maple('mtaylor',f,'[x=a]',3)
F =
(a^2-2*a)*exp(-a^2-y^2-a*y)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-2*a-y)+(2*a-2)*exp(-a^2-y^2-a*y))*(x-a)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-1+2*a^2+2*a*y+1/2*y^2)+exp(-a^2-y^2-a*y)+(2*a-2)*exp(-a^2-y^2-a*y)*(-2*a-y))*(x-a)^2
3.2.2 Fourier 級(jí)數(shù)展開
function [A,B,F]=fseries(f,x,n,a,b)
if nargin==3, a=-pi; b=pi; end
L=(b-a)/2;
if a+b, f=subs(f,x,x+L+a); end％變量區(qū)域互換
A=int(f,x,-L,L)/L; B=[]; F=A/2; %計(jì)算a0
for i=1:n
an=int(f*cos(i*pi*x/L),x,-L,L)/L;
bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=[A, an]; B=[B,bn];
F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L);
end
if a+b, F=subs(F,x,x-L-a); end ％換回變量區(qū)域
例:
>> syms x; f=x*(x-pi)*(x-2*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;? % 定義方波信號(hào)
>> xx=[-pi:pi/200:pi]; xx=xx(xx~=0); xx=sort([xx,-eps,eps]); % 剔除零點(diǎn)
>> yy=subs(f,x,xx); plot(xx,yy,'r-.'), hold on? % 繪制出理論值并保持坐標(biāo)系
>> 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)
……
3.2.3 級(jí)數(shù)求和的計(jì)算
是在符號(hào)工具箱中提供的
例:計(jì)算
>> format long; sum(2.^[0:63])?? %數(shù)值計(jì)算
ans =
1.844674407370955e+019
>> sum(sym(2).^[0:200]) % 或 syms k; symsum(2^k,0,200)
％把2定義為符號(hào)量可使計(jì)算更精確
ans =
3213876088517980551083924184682325205044405987565585670602751
>> syms k; symsum(2^k,0,200)
ans =
3213876088517980551083924184682325205044405987565585670602751
例:試求解無(wú)窮級(jí)數(shù)的和
>> syms n; s=symsum(1/((3*n-2)*(3*n+1)),n,1,inf)
%采用符號(hào)運(yùn)算工具箱
s =
1/3
>> m=1:10000000; s1=sum(1./((3*m-2).*(3*m+1)));%數(shù)值計(jì)算方法,雙精度有效位16,“大數(shù)吃小數(shù)”,無(wú)法精確
>> format long; s1 % 以長(zhǎng)型方式顯示得出的結(jié)果
s1 =
0.33333332222165
例:求解
>> syms n x
>> s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf);
>> simple(s1)? % 對(duì)結(jié)果進(jìn)行化簡(jiǎn)，MATLAB 6.5 及以前版本因本身 bug 化簡(jiǎn)很麻煩
ans =
log((((2*x+1)^2)^(1/2)+1)/(((2*x+1)^2)^(1/2)-1))
%實(shí)際應(yīng)為log((x+1)/x)
例:求
>> syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf)
ans =
eulergamma
>> vpa(ans, 70)? % 顯示 70 位有效數(shù)字
ans =
.5772156649015328606065120900824024310421593359399235988057672348848677
符號(hào)函數(shù)計(jì)算器
單變量符號(hào)函數(shù)計(jì)算器
Taylor 逼近計(jì)算器
單變量符號(hào)函數(shù)計(jì)算器(1/3)
在命令窗口中執(zhí)行 funtool 即可調(diào)出單變量符號(hào)函數(shù)計(jì)算器。單變量符號(hào)函數(shù)計(jì)算器用于對(duì)單變量函數(shù)進(jìn)行操作，可以對(duì)符號(hào)函數(shù)進(jìn)行化簡(jiǎn)、求導(dǎo)、繪制圖形等。該工具的界面如圖所示。
單變量符號(hào)函數(shù)計(jì)算器(2/3)
1．輸入框的功能
如圖：
單變量符號(hào)函數(shù)計(jì)算器(3/3)
單變量符號(hào)函數(shù)計(jì)算器應(yīng)用示例
Taylor 逼近計(jì)算器
Taylor 逼近計(jì)算器用于實(shí)現(xiàn)函數(shù)的 taylor 逼近。在命令窗口中輸入 taylortool，調(diào)出Taylor 逼近計(jì)算器，界面及功能如圖。
MAPLE 函數(shù)的調(diào)用
maple 函數(shù)的使用
mfun 函數(shù)的使用
maple 函數(shù)的使用
maple 是符號(hào)工具箱中的一個(gè)通用命令，使用它可以實(shí)現(xiàn)對(duì) MAPLE 中大部分函數(shù)的調(diào)用。其使用格式為：
r = maple('statement')，其中 statement 為符合 MAPLE 語(yǔ)法的可執(zhí)行語(yǔ)句的字符串，該命令將 statement 傳遞給 MAPLE，該命令的輸出結(jié)果也符合 MAPLE 的語(yǔ)法；
r = maple('function',arg1,arg2,...)，該函數(shù)調(diào)用引號(hào)中的函數(shù)，并接受指定的參數(shù)，相當(dāng)于 MAPLE 語(yǔ)句 function(arg1,arg2,...)；
[r, status] = maple(...)，返回函數(shù)的運(yùn)行狀態(tài)，如果函數(shù)運(yùn)行成功，則 status 為 0，r 為運(yùn)行結(jié)果；如果函數(shù)運(yùn)行失敗，則 status 為一個(gè)正數(shù)，r 為相應(yīng)的錯(cuò)誤信息；
maple('traceon') 或者 maple trace on，輸出 MAPLE? 函數(shù)運(yùn)行中的所有子表達(dá)式和運(yùn)行結(jié)果；
maple('traceoff') 或 maple trace off，不顯示中間過(guò)程。
mfun 函數(shù)的使用
mfun 函數(shù)用于對(duì) maple 函數(shù)進(jìn)行數(shù)字評(píng)估。該函數(shù)的調(diào)用格式為：
Y = mfun('function',par1,par2,par3,par4)。
該語(yǔ)句對(duì)指定的數(shù)學(xué)函數(shù)進(jìn)行評(píng)估。其中 'function' 指定待評(píng)估的函數(shù)，par1、par2 等為 'function' 的參數(shù)，為待評(píng)估的數(shù)值，其維數(shù)有 'function' 函數(shù)的參數(shù)類型確定。在該語(yǔ)句中最多可以設(shè)置四個(gè)參數(shù)，最后一個(gè)參數(shù)可以為矩陣。
用戶可以通過(guò) help mfunlist 查看 MATLAB 中 mfun 可以調(diào)用的函數(shù)列表，另外，可以通過(guò) mhelp function 查看指定函數(shù)的相關(guān)信息。
3.3.1 數(shù)值微分算法
向前差商公式：
向后差商公式
兩種中心公式：
3.3.2 中心差分方法及其 MATLAB 實(shí)現(xiàn)
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;
％數(shù)值計(jì)算diff(X)表示數(shù)組X相鄰兩數(shù)的差
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:end-L0); dx=([1:length(dy)]+L0-2-(n>2))*Dt;
調(diào)用格式：
y為 等距實(shí)測(cè)數(shù)據(jù)， dy為得出的導(dǎo)數(shù)向量， dx為相應(yīng)的自變量向量，dy、dx的數(shù)據(jù)比y短 。
例：
求導(dǎo)數(shù)的解析解,再用數(shù)值微分求取原函數(shù)的1～4 階導(dǎo)數(shù)，并和解析解比較精度。
>> h=0.05; x=0:h:pi;
>> syms x1; y=sin(x1)/(x1^2+4*x1+3);
% 求各階導(dǎo)數(shù)的解析解與對(duì)照數(shù)據(jù)
>> 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);?? % 生成已知數(shù)據(jù)點(diǎn)
>> [y1,dx1]=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,':');
>> [y2,dx2]=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,y2,':')
>> [y3,dx3]=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,':');
>> [y4,dx4]=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,':')
求最大相對(duì)誤差：
>> norm((y4-…
f4(4:60))./f4(4:60))
ans =
3.5025e-004
3.3.3? 用插值、擬合多項(xiàng)式的求導(dǎo)數(shù)
基本思想：當(dāng)已知函數(shù)在一些離散點(diǎn)上的函數(shù)值時(shí)，該函數(shù)可用插值或擬合多項(xiàng)式來(lái)近似，然后對(duì)多項(xiàng)式進(jìn)行微分求得導(dǎo)數(shù)。
選取x=0附近的少量點(diǎn)
進(jìn)行多項(xiàng)式擬合或插值
g(x)在x=0處的k階導(dǎo)數(shù)為
通過(guò)坐標(biāo)變換用上述方法計(jì)算任意x點(diǎn)處的導(dǎo)數(shù)值
令
將g(x)寫成z的表達(dá)式
導(dǎo)數(shù)為
可直接用??????? 擬合節(jié)點(diǎn)???????????????? 得到系數(shù)
d=polyfit(x-a,y,length(xd)-1)
例：數(shù)據(jù)集合如下：
xd: 0?????????? 0.2000? 0.4000? 0.6000? 0.8000? 1.000
yd: 0.3927? 0.5672? 0.6982? 0.7941? 0.8614? 0.9053
計(jì)算x=a=0.3處的各階導(dǎo)數(shù)。
>>? xd=[ 0? 0.2000? 0.4000? 0.6000? 0.8000? 1.000];
>>? yd=[0.3927? 0.5672? 0.6982? 0.7941? 0.8614? 0.9053];
>> a=0.3;L=length(xd);
>> d=polyfit(xd-a,yd,L-1);fact=[1];
>>? for k=1:L-1;fact=[factorial(k),fact];end
>>? deriv=d.*fact
deriv =
1.8750?? -1.3750??? 1.0406?? -0.9710??? 0.6533??? 0.6376
建立用擬合(插值)多項(xiàng)式計(jì)算各階導(dǎo)數(shù)的poly_drv.m
function der=poly_drv(xd,yd,a)
m=length(xd)-1;
d=polyfit(xd-a,yd,m);
c=d(m:-1:1);?? ％去掉常數(shù)項(xiàng)
fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);end
der=c.*fact;
例：
>>? xd=[ 0? 0.2000? 0.4000? 0.6000? 0.8000? 1.000];
>> yd=[0.3927? 0.5672? 0.6982? 0.7941? 0.8614? 0.9053];
>> a=0.3;?? der=poly_drv(xd,yd,a)
der =
0.6533?? -0.9710??? 1.0406?? -1.3750??? 1.8750
3.3.4 二元函數(shù)的梯度計(jì)算
格式：
若z矩陣是建立在等間距的形式生成的網(wǎng)格基礎(chǔ)上，則實(shí)際梯度為
例：
計(jì)算梯度，繪制引力線圖：
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z);
>> fx=fx/0.2; fy=fy/0.2;
>> contour(x,y,z,30);
>> hold on;
>> quiver(x,y,fx,fy)
%繪制等高線與
引力線圖
繪制誤差曲面:
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.08])
>> figure;? surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.11])
％建立一個(gè)新圖形窗口
為減少誤差，對(duì)網(wǎng)格加密一倍：
>> [x,y]=meshgrid(-3:.1:3,-2:.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z); fx=fx/0.1; fy=fy/0.1;
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.02])
>> figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.06])
3.4 數(shù)值積分問(wèn)題 4.3.1 由給定數(shù)據(jù)進(jìn)行梯形求積
格式： S=trapz(x,y)
例：
>> x1=[0:pi/30:pi]'; y=[sin(x1) cos(x1) sin(x1/2)];
>> x=[x1 x1 x1]; S=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
S =
1.9982??? 0.0000??? 1.9995
>> S1=trapz(x1,y)???? % 得出和上述完全一致的結(jié)果
S1 =
1.9982??? 0.0000??? 1.9995
例：
畫圖
>> x=[0:0.01:3*pi/2, 3*pi/2];? % 這樣賦值能確保 3*pi/2點(diǎn)被包含在內(nèi)
>> y=cos(15*x); plot(x,y)
% 求取理論值
>> syms x, A=int(cos…
(15*x),0,3*pi/2)
A =
1/15
隨著步距h的減小，計(jì)算精度逐漸增加：
>> h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; v=[];
>> for h=h0,
x=[0:h:3*pi/2, 3*pi/2]; y=cos(15*x); I=trapz(x,y);
v=[v; h, I, 1/15-I ];
end
>> format long; v
v =
0.10000000000000?? 0.05389175150076?? 0.01277491516591
0.01000000000000?? 0.06654169546584?? 0.00012497120083
0.00100000000000?? 0.06666541668004?? 0.00000124998663
0.00010000000000?? 0.06666665416667?? 0.00000001250000
0.00001000000000?? 0.06666666654167?? 0.00000000012500
0.00000100000000?? 0.06666666666542?? 0.00000000000125
3.4.2 單變量數(shù)值積分問(wèn)題求解
梯形公式
格式：(變步長(zhǎng))(Fun:函數(shù)的字符串變量)
y=quad(Fun,a,b)??? y=quadl(Fun,a,b) % 求定積分
y=quad(Fun,a,b,??? ) y=quadl(Fun,a,b,??? ) %限定精度的定積分求解，默認(rèn)精度為10－6。后面函數(shù)算法更精確，精度更高。
例：
函數(shù)定義被積函數(shù)：
>> y=quad('c3ffun',0,1.5)
y =
0.9661
用 inline 函數(shù)定義被積函數(shù)：
>> f=inline('2/sqrt(pi)*exp(-x.^2)','x');
>> y=quad(f,0,1.5)
y =
0.9661
運(yùn)用符號(hào)工具箱：
>> syms x, y0=vpa(int(2/sqrt(pi)*exp(-x^2),0,1.5),60)
y0 =
.966105146475310713936933729949905794996224943257461473285749
>> y=quad(f,0,1.5,1e-20)??? % 設(shè)置高精度，但該方法失效
提高求解精度：
>> y=quadl(f,0,1.5,1e-20)
y =
0.9661
>> abs(y-y0)
ans =
.6402522848913892e-16
>> format long? ％16位精度
>> y=quadl(f,0,1.5,1e-20)
y =
0.96610514647531
例：求解
繪制函數(shù):
>> x=[0:0.01:2, 2+eps:0.01:4,4];
>> y=exp(x.^2).*(x<=2)+80./(4-sin(16*pi*x)).*(x>2);
>> y(end)=0;
>> x=[eps, x];
>> y=[0,y];
>> fill(x,y,'g')
％為減少視覺(jué)上的誤
差，對(duì)端點(diǎn)與間斷點(diǎn)
(有跳躍)進(jìn)行處理。
調(diào)用quad( ):
>> f=inline('exp(x.^2).*(x<=2)+80*(x>2)./(4-sin(16*pi*x))','x');
>> I1=quad(f,0,4)
I1 =
57.76435412500863
調(diào)用quadl( ):
>> I2=quadl(f,0,4)
I2 =
57.76445016946768
>> syms x; I=vpa(int(exp(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4))
I =
57.764450125053010333315235385182
3.4.3?? Gauss求積公式
為使求積公式得到較高的代數(shù)精度
對(duì)求積區(qū)間[a,b]，通過(guò)變換
有
以n=2的高斯公式為例：
function g=gauss2(fun,a,b)
h=(b-a)/2;
c=(a+b)/2;
x=[h*(-0.7745967)+c, c, h*0.7745967+c];
g=h*(0.55555556*(gaussf(x(1))+gaussf(x(3)))+0.88888889*gaussf(x(2)));
function y=gaussf(x)
y=cos(x);
>> gauss2('gaussf',0,1)
ans =
0.8415
3.4.4 雙重積分問(wèn)題的數(shù)值解
矩形區(qū)域上的二重積分的數(shù)值計(jì)算
格式：
矩形區(qū)域的雙重積分：
y=dblquad(Fun,xm,xM,ym,yM)
限定精度的雙重積分：
y=dblquad(Fun,xm,xM,ym,yM，?? )
例：求解
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');
>> y=dblquad(f,-2,2,-1,1)
y =
1.57449318974494
任意區(qū)域上二元函數(shù)的數(shù)值積分 (調(diào)用工具箱NIT)，該函數(shù)指定順序先x后y.
解析解方法：
>> syms x y
>> i1=int(exp(-x^2/2)*sin(x^2+y), y, -sqrt(1-x^2/2), sqrt(1-x^2/2));
>> int(i1, x, -1/2, 1)
Warning: Explicit integral could not be found.
> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58
ans =
int(2*exp(-1/2*x^2)*sin(x^2)*sin(1/2*(4-2*x^2)^(1/2)), x = -1/2 .. 1)
>> vpa(ans)
ans =
.41192954617629511965175994017601
例：計(jì)算單位圓域上的積分：
先把二重積分轉(zhuǎn)化：
對(duì)x是不可積的，故調(diào)用解析解方法不會(huì)得出結(jié)果，而數(shù)值解求解不受此影響。
>> fh=inline('sqrt(1-y.^2)','y');? % 內(nèi)積分上限
>> fl=inline('-sqrt(1-y.^2)','y'); % 內(nèi)積分下限
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');? %交換順序的被積函數(shù)
>> I=quad2dggen(f,fl,fh,-1,1,eps)
Integral did not converge--singularity likely
I =
0.53686038269795
3.4.5 三重定積分的數(shù)值求解
格式:
I=triplequad(Fun,xm,xM,ym,yM, zm,zM，??? ,@quadl)
其中@quadl為具體求解一元積分的數(shù)值函數(shù),也可選用@quad或自編積分函數(shù),但調(diào)用格式要與quadl一致。
例：
>> triplequad(inline('4*x.*z.*exp(-x.*x.*y-z.*z)', …
'x','y','z'), 0, 1, 0, pi, 0, pi,1e-7,@quadl)
ans =
1.7328
3.5 曲線積分與曲面積分的計(jì)算
3.5.1 曲線積分及MATLAB求解
第一類曲線積分
起源于對(duì)不均勻分布的空間曲線總質(zhì)量的求取.設(shè)空間曲線L的密度函數(shù)為f(x,y,z),則其總質(zhì)量
其中s為曲線上某點(diǎn)的弧長(zhǎng),又稱這類曲線積分為對(duì)弧長(zhǎng)的曲線積分.
數(shù)學(xué)表示
若
弧長(zhǎng)表示為
例：
>> syms t; syms a positive; x=a*cos(t); y=a*sin(t); z=a*t;
>> I=int(z^2/(x^2+y^2)*sqrt(diff(x,t)^2+diff(y,t)^2+ diff(z,t)^2),t,0,2*pi)
I =
8/3*pi^3*a*2^(1/2)
>> pretty(I)
3???? 1/2
8/3 pi? a 2
例：
>> x=0:.001:1.2; y1=x; y2=x.^2; plot(x,y1,x,y2)
％繪出兩條曲線
>> syms x; y1=x; y2=x^2; I1=int((x^2+y2^2)*sqrt(1+diff(y2,x)^2),x,0,1);
>> I2=int((x^2+y1^2)*sqrt(1+diff(y1,x)^2),x,1,0); I=I2+I1
I =
-2/3*2^(1/2)+349/768*5^(1/2)+7/512*log(-2+5^(1/2))
3.5.1.2 第二類曲線積分
又稱對(duì)坐標(biāo)的曲線積分，起源于變力
沿曲線??? 移動(dòng)時(shí)作功的研究
曲線????? 亦為向量,若曲線可以由參數(shù)方程表示
則兩個(gè)向量的點(diǎn)乘可由這兩個(gè)向量直接得出.
例：求曲線積分
>> syms t; syms a positive; x=a*cos(t); y=a*sin(t);
>> F=[(x+y)/(x^2+y^2),-(x-y)/(x^2+y^2)]; ds=[diff(x,t);diff(y,t)];
>> I=int(F*ds,t,2*pi,0)? % 正向圓周
I =
2*pi
例：
>> syms x; y=x^2; F=[x^2-2*x*y,y^2-2*x*y]; ds=[1; diff(y,x)];
>> I=int(F*ds,x,-1,1)
I =
-14/15
3.5.2曲面積分與MATLAB語(yǔ)言求解3.5.2.1 第一類曲面積分
其中???? 為小區(qū)域的面積，故又稱為對(duì)面積的曲面積分。曲面??? 由?????????????? 給出，則該積分可轉(zhuǎn)換成x-y平面的二重積分為
例：
％四個(gè)平面，其中三個(gè)被積函數(shù)的值為0，只須計(jì)算一個(gè)即可。
>> syms x y; syms a positive; z=a-x-y;
>> I=int(int(x*y*z*sqrt(1+diff(z,x)^2+ diff(z,y)^2),y,0,a-x),x,0,a)
I =
1/120*3^(1/2)*a^5
若曲面由參數(shù)方程
曲面積分
例：
>> syms u v; syms a positive;
>> x=u*cos(v); y=u*sin(v); z=v;f=x^2*y+z*y^2;
>> E=simple(diff(x,u)^2+diff(y,u)^2+diff(z,u)^2);
>> F=diff(x,u)*diff(x,v)+diff(y,u)*diff(y,v)+diff(z,u)* diff(z,v);
>> G=simple(diff(x,v)^2+diff(y,v)^2+diff(z,v)^2);
>> I=int(int(f*sqrt(E*G-F^2),u,0,a),v,0,2*pi)
I =
1/4*a*(a^2+1)^(3/2)*pi^2+1/8*log(-a+(a^2+1)^(1/2)) *pi^2-1/8*(a^2+1)^(1/2)*a*pi^2
3.5.2.2 第二類曲面積分
又稱對(duì)坐標(biāo)的曲面積分
可轉(zhuǎn)化成第一類曲面積分
若曲面由參數(shù)方程給出
例：
的上半部，且積分沿橢球面的上面。
％引入?yún)?shù)方程?? x=a*sin(u)*cos(v); y=b*sin(u)*sin(v);? z=c*cos(u)， u[0,pi/2], v[0,2*pi].
>> syms u v; syms a b c positive;
>> x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u);
>> A=diff(y,u)*diff(z,v)-diff(z,u)*diff(y,v);
>> I=int(int(x^3*A,u,0,pi/2),v,0,2*pi)
I =
2/5*pi*a^3*c*b
相關(guān)PPT
《中科院的matlab課件》是由用戶諍友于2017-08-20上傳，屬于課件PPT。
總結(jié)
以上是生活随笔為你收集整理的中科院的matlab课件,中科院的matlab课件的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： php客户反馈发邮箱,PHP实现通过ge
- 下一篇： matlab 折线图_MATLAB作图实