Question

有关于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].

麻烦写出详细的程序步骤 谢谢

Answer 1

第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''')

 

Answer 2

这应该是第一类曲线和曲面积分，
1
积分曲线是个椭圆，可以设出参数方程x=cost, y=4sint

积分m=∫(0->2π){2[cost-(4sint)^2]^3√[(-sint)^2+(4cost)^2]} dt

2
积分曲面上，z'x=6(x-y^2)^2, z'y=6(x-y^2)^2*(-2y)=-12y(x-y^2)^2
积分m=∫∫{x^2+y^2+[2(x-y^2)^3]^2}√[1+(z'x)^2+(z'y)^2] dxdy
带入z'x和z'y即可。。

写程序不会，不过写到这一部，带入wolframalpha也能直接算出结果了。。