有关于Matlab的两道计算题1.某个函数（类似于z=2(x-y^2)^3）在x^2+(y/4)^2=1这个平面曲线上的积分.2.f=x^2+y^2+z^2在函数z=2(x-y^2)^3的空间曲面上的积分.其中x跟y的积分区间均为[-6,6].

有关于Matlab的两道计算题1.某个函数（类似于z=2(x-y^2)^3）在x^2+(y/4)^2=1这个平面曲线上的积分.2.f=x^2+y^2+z^2在函数z=2(x-y^2)^3的空间曲面上的积分.其中x跟y的积分区间均为[-6,6].
有关于Matlab的两道计算题
1.某个函数（类似于z=2(x-y^2)^3）在x^2+(y/4)^2=1这个平面曲线上的积分.2.f=x^2+y^2+z^2在函数z=2(x-y^2)^3的空间曲面上的积分.其中x跟y的积分区间均为[-6,6].

有关于Matlab的两道计算题1.某个函数（类似于z=2(x-y^2)^3）在x^2+(y/4)^2=1这个平面曲线上的积分.2.f=x^2+y^2+z^2在函数z=2(x-y^2)^3的空间曲面上的积分.其中x跟y的积分区间均为[-6,6].
第1题
属于第一类曲线积分,计算原理可以参考附件ppt的92-93页（内容源自薛定宇《高等数学问题的MATLAB求解》,下同）.
使用MATLAB符号积分函数int计算,但无法求出解析解,可通过vpa得到所需精度的数值结果.
代码如下：
syms t
x = cos(t);
y = 4*sin(t);
z = 2*(x-y^2)^3;
F = z * sqrt(diff(x,t)^2+diff(y,t)^2);
I = int(F,t,0,2*pi);
vpa(I)
得到结果为
ans = 
-25865.95635790842439629780459932
也可以直接使用 数值方法求解,把上述代码最后两行改成
quadl(inline(vectorize(F)),0,2*pi)
或
quadl(matlabFunction(F),0,2*pi)
得到结果为
ans =
 -2.5866e+004
 
第2题
属于第一类曲面积分,计算原理可以参考附件ppt的98-99页.
用MATLAB计算可以使用符号积分函数int,但无法求出解析解,然后使用vpa得到所需精度的数值结果.
代码如下：
syms x y
z = 2*(x-y^2)^3;
f = x^2+y^2+z^2;
F = f * sqrt(1+diff(z,x)^2+diff(z,y)^2);
I = int(int(F,x,-6,6),y,-6,6)
vpa(I)
得到结果：
ans = 
8829165547218478.3258650115855033
 
需要说明的是,使用int+vpa的方法耗时较长,而且不同MATLAB版本上的情况可能存在差别.我在2008b上计算两重int耗时超过半小时,得到最后结果共耗时秒（和具体硬件、软件环境相关,仅供参考）.
 
如果采用数值积分,由于在x=-6,y=±6附近函数的值较大（下图为把曲面积分转换成普通二重积分的被积函数曲面图形,使用 ezmesh(F,[-6 6]) 绘制）,如果用默认的误差限（1e-6）,计算时会出现大量的“Maximum function count exceeded; singularity likely”警告信息,计算精度无法得到保证.

使用大一些的误差限能避免该问题：
>> dblquad(inline(vectorize(F)),-6,6,-6,6,0.1,@quadl)
ans =
  8.8292e+015
或者,使用匿名函数的形式比起inline函数速度更快一些：
dblquad(matlabFunction(F),-6,6,-6,6,0.1,@quadl)
 
第3题（a）
属于常微分方程的初值问题,可用ode系列函数求解.
微分方程看不太清楚,等号右边第1项是不是y的二阶导数?
dy = inline('[y(2); y(3); 2*y(3)+6*x*y(1)]','x','y');
y0 = [1;1;1];
[T1, Y1] = ode45(dy,[0 10], y0);
[T2, Y2] = ode45(dy,[0 -10], y0);
T = [T2(end:-1:2); T1];
Y = [Y2(end:-1:2, :); Y1];
plot(T,Y)
得到结果看上去是发散的（不确定方程是否有问题）：

 
第3题（b）
属于边值问题（BVP）,可用bvp4c求
dy = inline('[y(2); y(3); 2*y(3)+6*x*y(1)]','x','y');
bc = inline('[ya(1)-2; yb(1); yb(2)]','ya','yb');
solinit = bvpinit(linspace(0,5,5),[0 0 0]);
sol = bvp4c(dy,bc,solinit);
x = linspace(0,5);
y = deval(sol,x);
subplot 211
plot(x,y(1,:));
ylabel('y')
subplot 212
plot(x,y(2,:));
xlabel('x')
ylabel('y''')