Win8+Matlab7.1测试例子（20150109、0116、0124、0202、0207、0217）

20150403-20150412补充：

用到的函数文件mo.m
function rchamo=mo(r) %求向量r的模
sum=0;
for i=1:3
sum=sum+r(i)^2;
end
rchamo=simple(sum^(1/2));

求曲线在任一点处的曲率与挠率
空间曲线曲率与挠率的计算：
求曲线r=(t,t^2/2,t^3/3)在点t=1上的曲率向量k和曲率|k|。
>> syms t;r=[t t^2/2 t^3/3];r1=diff(r)
 
r1 =
 
[   1,   t, t^2]
 
 
>> r1mo=mo(r1)
 
r1mo =
 
(1+t^2+t^4)^(1/2)
 
 
>> T=r1/r1mo
 
T =
 
[   1/(1+t^2+t^4)^(1/2),   t/(1+t^2+t^4)^(1/2), t^2/(1+t^2+t^4)^(1/2)]
 
 
>> T1=diff(T)
 
T1 =
 
[                          -1/2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3),     1/(1+t^2+t^4)^(1/2)-1/2*t/(1+t^2+t^4)^(3/2)*(2*t+4*t^3), 2*t/(1+t^2+t^4)^(1/2)-1/2*t^2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3)]
 
 
>> k=T1/r1mo
 
k =
 
[                                                  -1/2/(1+t^2+t^4)^2*(2*t+4*t^3),     (1/(1+t^2+t^4)^(1/2)-1/2*t/(1+t^2+t^4)^(3/2)*(2*t+4*t^3))/(1+t^2+t^4)^(1/2), (2*t/(1+t^2+t^4)^(1/2)-1/2*t^2/(1+t^2+t^4)^(3/2)*(2*t+4*t^3))/(1+t^2+t^4)^(1/2)]
 
 
>> t=1;kzhi=eval(k)

kzhi =

   -0.3333    0.0000    0.3333 %曲率向量k=(-1/3,0,1/3)
>> kzhimo=sqrt(kzhi*kzhi')

kzhimo =

    0.4714 %曲率|k|=sqrt(2)/3
求曲线r=(3t-t^3,3t^2,3t+t^3)的曲率和挠率。
>> syms t;r=[3*t-t^3 3*t^2 3*t+t^3];r1=diff(r)
 
r1 =
 
[ 3-3*t^2,     6*t, 3+3*t^2]
 
 
>> r2=diff(r1)
 
r2 =
 
[ -6*t,    6,  6*t]
 
 
>> r1mo=mo(r1)
 
r1mo =
 
3*2^(1/2)*(1+t^2)
 
 
>> r1r2=cross(r1,r2)
 
r1r2 =
 
[                    18*t^2-18, -6*(3+3*t^2)*t-6*(3-3*t^2)*t,                    18+18*t^2]
 
 
>> r1r2mo=mo(r1r2)
 
r1r2mo =
 
18*2^(1/2)*(1+t^2)
 
 
>> k=r1r2mo/r1mo^3
 
k =
 
1/3/(1+t^2)^2
 
 
>> r3=diff(r2)
 
r3 =
 
[ -6,  0,  6]
 
 
>> r123=dot(cross(r1,r2),r3)
 
r123 =
 
216
 
 
>> tau=r123/r1r2mo^2
 
tau =
 
1/3/(1+t^2)^2 %注意，沿此曲线，|k|=tau=1/(3*(1+t^2)^2)

用符号积分求解二重积分∫(-1,2)∫(y*y,y+2)xydxdy
>> syms x y;%定义两个符号变量
>> a=int(int(x*y,x,y*y,y+2),y,-1,2);%积分限x：y*y,y+2   ,y:,-1,2
>> b=simple(a);%化简
>> c=vpa(b,4)%得到4位近似解，也可以任意N位解
 
c =
 
5.625
a =
 
45/8
b =
 
45/8
>> int(int(x*y,x,y*y,y+2),y,-1,2)
 
ans =
 
45/8
dblquad计算二重数值积分
>> f=@(x,y)x.*y;
>> dblquad(f,y*y,y+2,-1,2)
出错
>> dblquad(f,pi/4,1,2,4)

ans =

    1.1494

多重积分的MATLAB实现
http://wenku.baidu.com/link?url=hk9pQjn2em9epQpLXaKreqU3FwFKI7SB1g5-NE8klWTE2ywsOIY2u8s4pDy4w8LplEg_EFgtuza2dwZEu9_sGZn4bm59UVljxEkDSrEC-6m

将积分分为定积分、二重积分、三重积分、第一型曲线积分、第二型曲线积分、第一型曲面积分和第二型曲面积分等七种积分
例1计算二重积分，积分对象是(x^2+y^2-x)dxdy=(x^2+y^2-x)dx∧dy，积分区域是由直线y=2,y=x,y=2x所围成的闭区域
>> int(int(x^2+y^2-x,x,y/2,y),y,0,2)
 
ans =
 
13/6
例9计算三重积分，积分对象是xyzdxdydz=xyzdx∧dy∧dz，积分区域是球面x^2+y^2+z^2=1及三个坐标面所围成的在第一象限内的区域
>> syms x y z;int(int(int(x*y*z,z,0,sqrt(1-x^2-y^2)),y,0,sqrt(1-x^2)),x,0,1)
 
ans =
 
1/48
计算定积分，函数y=x^2在(1,2)区间内的投影面积
>> f=@(x)x.^2;quad(f,1,2) %算子（泛函）quad(f,a,b)表示在区间(a,b)内对x求积分

ans =

    2.3333
MATLAB积分
http://wenku.baidu.com/link?url=vyGVVVoPFxXtd-CLC7dwdqlX6S140Q0gAhDNnUwU3i6Ilm57aR7PnKrr5Ja6uiymFInaof6h0hzA1vhyiZCJGVx1nYzFqBFwGANTTz_gX7K
triplequad计算三重数值积分
>> f=@(x,y,z)(y*sin(x)+z*cos(x));triplequad(f,0,pi,0,1,-1,1)

ans =

    2.0000
实验七+曲线、曲面积分及其应用
http://wenku.baidu.com/link?url=NZJ54kMV82K117znxE62kN5UprabK758Xj0GfkH--NKU_aTeVYrLl-TYKIy0xXUIUQ_vWCaboouputuLfiRE7dPMh1S3FqG8HF48G9ufzlW
用第二型曲线积分求椭圆的面积（利用格林公式计算平面图形的面积）
>> syms x y t dx dy a b
>> x=a*cos(t);y=b*sin(t);
>> dx=diff(x);dy=diff(y);
>> int(x*dy-y*dx,t,0,2*pi)
 
ans =
 
2*a*b*pi
 
 
>> (1/2)*int(x*dy-y*dx,t,0,2*pi)
 
ans =
 
a*b*pi
——椭圆x^2/a^2+y^2/b^2=1的面积为abpi

用第一型曲线积分求椭圆的周长
>> ds=sqrt(dx^2+dy^2)
 
ds =
 
(a^2*sin(t)^2+b^2*cos(t)^2)^(1/2)
 
 
>> int(ds,t,0,2*pi)
 
ans =
 
4*EllipticE((-(-a^2+b^2)/a^2)^(1/2))*(b^2/a^2)^(1/2)*a^2/(b^2)^(1/2)
>> simple(ans)

 ans =
 
4*EllipticE(1/a*(a^2-b^2)^(1/2))*a

雅可比用q来表示模和周期，例如
k=4sqrt(q)[(1+q^2)(1+q^4)(1+q^6)…/((1+q)(1+q^3)(1+q^5)…)]^4。
k=kq(0.207879576350762)=0.985171431009416
>> syms kq1 q;kq1=maple('product(((1+q^(2*n))/(1+q^(2*n-1))),n=1..100)');q=0.207879576350762;kq=eval(kq1);k=4*sqrt(q)*kq^4

k =

    0.9852

Matlab在数论中的应用
>> ISPRIME(3)

ans =

     1

>> isprime(3)

ans =

     1

>> ISPRIME(30)

ans =

     0

>> isprime(30)

ans =

     0
isprime测试n是否为素数

>> maple('gcd(14, 21)')

ans =

7
>> maple('ithprime(3)')

ans =

5
ithprime计算第n个素数
>> maple('prevprime(100)')

ans =

97
prevprime计算比n小的最大的素数
>> maple('ifactor(100)')

ans =

``(2)^2*``(5)^2
ifactor求整数n的素幂分解
>> maple('ifactors(100)')

ans =

[1, [[2, 2], [5, 2]]]
ifactors求整数n的所有素因子
>> maple('numtheory[phi](5)')

ans =

4
numtheory[phi]计算欧拉phi函数在n处的值
>> maple('numtheory[fermat](5)')

ans =

4294967297
numtheory[fermat]计算第n个费马数
>> maple('numtheory[invphi](12)')

ans =

[13, 21, 26, 28, 36, 42]
numtheory[invphi]计算使得phi(m)=n的正整数m
>> maple('numtheory[tau](12)')

ans =

6
numtheory[tau]计算n的所有正因子的个数
>> maple('numtheory[bigomega](12)')

ans =

3
numtheory[bigomega]计算n的素因子个数omega(n)的值
>> maple('numtheory[primroot](12)')

ans =

FAIL

>> maple('numtheory[primroot](7)')

ans =

3
numtheory[primroot]计算模n的最小原根
>> maple('numtheory[legendre](5,7)')

ans =

-1
numtheory[legendre]计算勒让德符号的值
>> maple('numtheory[jacobi](5,7)')

ans =

-1
numtheory[jacobi]计算雅可比符号的值
>> maple('with(GaussInt);')

ans =

[GIbasis, GIchrem, GIdivisor, GIfacpoly, GIfacset, GIfactor, GIfactors, GIgcd, GIgcdex, GIhermite, GIissqr,

GIlcm, GImcmbine, GInearest, GInodiv, GInorm, GInormal, GIorder, GIphi, GIprime, GIquadres, GIquo, GIrem,

GIroots, GIsieve, GIsmith, GIsqrfree, GIsqrt, GIunitnormal]
>> maple('GaussInt[GIprime](5+7*i);')

ans =

false

>> maple('GaussInt[GIprime](5+7*i)')

ans =

false
GaussInt[GIprime]判断高斯整数m是否为高斯素数

纪岗：Matlab语言在初等数论中的应用
http://www.docin.com/p-970648613.html
1.1求两个自然数的最大公约数的Matlab程序
a=input('请输入正整数a=');
b=input('请输入正整数b=');
yue=1;
for i=1:1:a
c=mod(a,i);
d=mod(b,i);
if c==0 && d==0
yue=i;
end
end
disp(['正整数',num2str(a),'与',num2str(b),'的最大公约数是',num2str(yue)])
请输入正整数a=18
请输入正整数b=24
正整数18与24的最大公约数是6
 >> bei=a*b/yue

bei =

    72
2.1生成小于n的所有素数表的Matlab程序——应该是小于等于
n=input('请输入正整数n=');
prime=[2 3 5];
for i=6:n
p=1;
for j=2:floor(sqrt(i))
if mod(i,j)==0
p=0;
break;
end
end
if p
prime=[prime i];%生成素数表
end
end
disp(['小于正整数',num2str(n),'的素数有',num2str(prime)])
请输入正整数n=30
小于正整数30的素数有2   3   5   7  11  13  17  19  23  29
请输入正整数n=30
小于正整数30的素数有2   3   5   7  11  13  17  19  23  29
>> find(prime==n)

ans =

   Empty mat