收藏 分销(赏)

数学实验课后习题解答.doc

上传人:xrp****65 文档编号:9434869 上传时间:2025-03-26 格式:DOC 页数:75 大小:1.31MB 下载积分:10 金币
下载 相关 举报
数学实验课后习题解答.doc_第1页
第1页 / 共75页
数学实验课后习题解答.doc_第2页
第2页 / 共75页


点击查看更多>>
资源描述
数学实验课后习题解答 配套教材:王向东 戎海武 文翰 编著数学实验 王汝军编写 实验一 曲线绘图 【练习与思考】 画出下列常见曲线的图形。 以直角坐标方程表示的曲线: 1. 立方曲线 clear; x=-2:0.1:2; y=x.^3; plot(x,y) 2. 立方抛物线 clear; y=-2:0.1:2; x=y.^3; plot(x,y) grid on 3. 高斯曲线 clear; x=-3:0.1:3; y=exp(-x.^2); plot(x,y); grid on %axis equal 以参数方程表示的曲线 4. 奈尔抛物线 clear; t=-3:0.05:3; x=t.^3;y=t.^2; plot(x,y) axis equal grid on 5. 半立方抛物线 clear; t=-3:0.05:3; x=t.^2;y=t.^3; plot(x,y) %axis equal grid on 6. 迪卡尔曲线 clear; a=3;t=-6:0.1:6; x=3*a*t./(1+t.^2); y=3*a*t.^2./(1+t.^2); plot(x,y) 7. 蔓叶线 clear; a=3;t=-6:0.1:6; x=3*a*t.^2./(1+t.^2); y=3*a*t.^3./(1+t.^2); plot(x,y) 8. 摆线 clear;clc; a=1;b=1; t=0:pi/50:6*pi; x=a*(t-sin(t)); y=b*(1-cos(t)); plot(x,y); axis equal grid on 9. 内摆线(星形线) clear; a=1; t=0:pi/50:2*pi; x=a*cos(t).^3; y=a*sin(t).^3; plot(x,y) 10. 圆的渐伸线(渐开线) clear; a=1; t=0:pi/50:6*pi; x=a*(cos(t)+t.*sin(t)); y=a*(sin(t)+t.*cos(t)); plot(x,y) grid on 11. 空间螺线 clear a=3;b=2;c=1; t=0:pi/50:6*pi; x=a*cos(t); y=b*sin(t); z=c*t; plot3(x,y,z) grid on 以极坐标方程表示的曲线: 12. 阿基米德线 clear; a=1; phy=0:pi/50:6*pi; rho=a*phy; polar(phy,rho,'r-*') 13. 对数螺线 clear; a=0.1; phy=0:pi/50:6*pi; rho=exp(a*phy); polar(phy,rho) 14. 双纽线 clear; a=1; phy=-pi/4:pi/50:pi/4; rho=a*sqrt(cos(2*phy)); polar(phy,rho) hold on polar(phy,-rho) 15. 双纽线 clear; a=1; phy=0:pi/50:pi/2; rho=a*sqrt(sin(2*phy)); polar(phy,rho) hold on polar(phy,-rho) 16. 四叶玫瑰线 clear;close a=1; phy=0:pi/50:2*pi; rho=a*sin(2*phy); polar(phy,rho) 17. 三叶玫瑰线 clear;close a=1; phy=0:pi/50:2*pi; rho=a*sin(3*phy); polar(phy,rho) 18. 三叶玫瑰线 clear;close a=1; phy=0:pi/50:2*pi; rho=a*cos(3*phy); polar(phy,rho) 实验二 极限与导数 【练习与思考】 1. 求下列各极限 (1) (2) (3) clear; syms n y1=limit((1-1/n)^n,n,inf) y2=limit((n^3+3^n)^(1/n),n,inf) y3=limit(sqrt(n+2)-2*sqrt(n+1)+sqrt(n),n,inf) y1 =1/exp(1) y2 =3 y3 =0 (4) (5) (6) clear; syms x ; y4=limit(2/(x^2-1)-1/(x-1),x,1) y5=limit(x*cot(2*x),x,0) y6=limit(sqrt(x^2+3*x)-x,x,inf) y4 =-1/2 y5 =1/2 y6 =3/2 (7) (8) (9) clear; syms x m y7=limit(cos(m/x),x,inf) y8=limit(1/x-1/(exp(x)-1),x,1) y9=limit(((1+x)^(1/3)-1)/x,x,0) y7 =1 y8 =(exp(1) - 2)/(exp(1) - 1) y9 =1/3 2. 考虑函数 作出图形,并说出大致单调区间;使用diff求,并求确切的单调区间。 clear;close; syms x; f=3*x^2*sin(x^3); ezplot(f,[-2,2]) grid on 大致的单调增区间:[-2,-1.7],[-1.3,1.2],[1.7,2]; 大致的单点减区间:[-1.7,-1.3],[1.2,1.7]; f1=diff(f,x,1) ezplot(f1,[-2,2]) line([-5,5],[0,0]) grid on axis([-2.1,2.1,-60,120]) f1 = 6*x*sin(x^3) + 9*x^4*cos(x^3) 用fzero函数找的零点,即原函数的驻点 x1=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-2,-1.7]) x2=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-1.7,-1.5]) x3=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[-1.5,-1.1]) x4=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',0) x5=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1,1.5]) x6=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1.5,1.7]) x7=fzero('6*x*sin(x^3) + 9*x^4*cos(x^3)',[1.7,2]) x1 = -1.9948 x2 = -1.6926 x3 = -1.2401 x4 = 0 x5 = 1.2401 x6 = 1.6926 x7 = 1.9948 确切的单调增区间:[-1.9948,-1.6926],[-1.2401,1.2401],[1.6926,1.9948] 确切的单调减区间:[-2,-1.9948],[-1.6926,-1.2401],[1.2401,1.6926],[1.9948,2] 3. 对于下列函数完成下列工作,并写出总结报告,评论极值与导数的关系, (i) 作出图形,观测所有的局部极大、局部极小和全局最大、全局最小值点的粗略位置; (iI) 求所有零点(即的驻点); (iii) 求出驻点处的二阶导数值; (iv) 用fmin求各极值点的确切位置; (v) 局部极值点与有何关系? (1) (2) (3) clear;close; syms x; f=x^2*sin(x^2-x-2) ezplot(f,[-2,2]) grid on f = x^2*sin(x^2 - x - 2) 局部极大值点为:-1.6,局部极小值点为为:-0.75,-1.6 全局最大值点为为:-1.6,全局最小值点为:-3 f1=diff(f,x,1) ezplot(f1,[-2,2]) line([-5,5],[0,0]) grid on axis([-2.1,2.1,-6,20]) f1 = 2*x*sin(x^2 - x - 2) + x^2*cos(x^2 - x - 2)*(2*x - 1) 用fzero函数找的零点,即原函数的驻点 x1=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-2,-1.2]) x2=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-1.2,-0.5]) x3=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[-0.5,1.2]) x4=fzero('2*x*sin(x^2-x-2)+x^2*cos(x^2-x-2)*(2*x-1)',[1.2,2]) x1 = -1.5326 x2 = -0.7315 x3 = -3.2754e-027 x4 = 1.5951 ff=@(x) x.^2.*sin(x.^2-x-2) ff(-2),ff(x1),ff(x2),ff(x3),ff(x4),ff(2) ff = @(x)x.^2.*sin(x.^2-x-2) ans = -3.0272 ans = 2.2364 ans = -0.3582 ans = -9.7549e-054 ans = -2.2080 ans = 0 实验三 级数 【练习与思考】 1. 用taylor命令观测函数的Maclaurin展开式的前几项, 然后在同一坐标系里作出函数和它的Taylor展开式的前几项构成的多项式函数的图形,观测这些多项式函数的图形向的图形的逼近的情况 (1) clear; syms x y=asin(x); y1=taylor(y,0,1) y2=taylor(y,0,5) y3=taylor(y,0,10) y4=taylor(y,0,15) x=-1:0.1:1; y=subs(y,x); y1=subs(y1,x); y2=subs(y2,x); y3=subs(y3,x); y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = 0 y2 = x^3/6 + x y3 = (35*x^9)/1152 + (5*x^7)/112 + (3*x^5)/40 + x^3/6 + x y4 = (231*x^13)/13312 + (63*x^11)/2816 + (35*x^9)/1152 + (5*x^7)/112 + (3*x^5)/40 + x^3/6 + x (2) clear; syms x y=atan(x);y1=taylor(y,0,3) y2=taylor(y,0,5),y3=taylor(y,0,10),y4=taylor(y,0,15) x=-1:0.1:1; y=subs(y,x);y1=subs(y1,x);y2=subs(y2,x); y3=subs(y3,x);y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = x y2 = x - x^3/3 y3 = x^9/9 - x^7/7 + x^5/5 - x^3/3 + x y4 = x^13/13 - x^11/11 + x^9/9 - x^7/7 + x^5/5 - x^3/3 + x (3) clear; syms x y=exp(x^2); y1=taylor(y,0,3) y2=taylor(y,0,5) y3=taylor(y,0,10) y4=taylor(y,0,15) x=-1:0.1:1; y=subs(y,x); y1=subs(y1,x); y2=subs(y2,x); y3=subs(y3,x); y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = x^2 + 1 y2 = x^4/2 + x^2 + 1 y3 = x^8/24 + x^6/6 + x^4/2 + x^2 + 1 y4 = x^14/5040 + x^12/720 + x^10/120 + x^8/24 + x^6/6 + x^4/2 + x^2 + 1 (4) clear; syms x y=sin(x)^2; y1=taylor(y,0,1) y2=taylor(y,0,5) y3=taylor(y,0,10) y4=taylor(y,0,15) x=-pi:0.1:pi; y=subs(y,x); y1=subs(y1,x); y2=subs(y2,x); y3=subs(y3,x); y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = 0 y2 = x^2 - x^4/3 y3 = - x^8/315 + (2*x^6)/45 - x^4/3 + x^2 y4 = (4*x^14)/42567525 - (2*x^12)/467775 + (2*x^10)/14175 - x^8/315 + (2*x^6)/45 - x^4/3 + x^2 (5) clear; syms x y=exp(x)/(1-x); y1=taylor(y,0,3) y2=taylor(y,0,5) y3=taylor(y,0,10) y4=taylor(y,0,15) x=-1:0.1:0; y=subs(y,x); y1=subs(y1,x); y2=subs(y2,x); y3=subs(y3,x); y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = (5*x^2)/2 + 2*x + 1 y2 = (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1 y3 = (98641*x^9)/36288 + (109601*x^8)/40320 + (685*x^7)/252 + (1957*x^6)/720 + (163*x^5)/60 + (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1 y4 = (47395032961*x^14)/17435658240 + (8463398743*x^13)/3113510400 + (260412269*x^12)/95800320 + (13563139*x^11)/4989600 + (9864101*x^10)/3628800 + (98641*x^9)/36288 + (109601*x^8)/40320 + (685*x^7)/252 + (1957*x^6)/720 + (163*x^5)/60 + (65*x^4)/24 + (8*x^3)/3 + (5*x^2)/2 + 2*x + 1 (6) clear; syms x y=log(x+sqrt(1+x^2)); y1=taylor(y,0,3) y2=taylor(y,0,5) y3=taylor(y,0,10) y4=taylor(y,0,15) x=-1:0.1:1; y=subs(y,x); y1=subs(y1,x); y2=subs(y2,x); y3=subs(y3,x); y4=subs(y4,x); plot(x,y,x,y1,':',x,y2,'-.',x,y3,'--',x,y4,':','linewidth',3) y1 = x y2 = x - x^3/6 y3 = (35*x^9)/1152 - (5*x^7)/112 + (3*x^5)/40 - x^3/6 + x y4 = (231*x^13)/13312 - (63*x^11)/2816 + (35*x^9)/1152 - (5*x^7)/112 + (3*x^5)/40 - x^3/6 + x 2. 求公式中的数的值. k=[4 5 6 7 8]; syms n symsum(1./n.^(2*k),1,inf) ans = [ pi^8/9450, pi^10/93555, (691*pi^12)/638512875, (2*pi^14)/18243225, (3617*pi^16)/325641566250] 3. 利用公式来计算的近似值。精确到小数点后100位,这时应计算到这个无穷级数的前多少项?请说明你的理由. 解:Matlab代码为 clear;clc;close epsl=1.0e-100; ep=1;fn=1;a=1;n=1; while ep>epsl a=a+fn; n=n+1; fn=fn/n; ep=fn; end fn vpa(a,100) n fn = 8.3482e-101 ans = 2.71828182845904553488480814849026501178741455078125 n = 70 精确到小数点后100位,这时应计算到这个无穷级数的前71项,理由是误差小于10的负100次方,需要最后一项小于10的负100次方,由上述循环知n=70时最后一项小于10的负100次方,故应计算到这个无穷级数的前71项. 4. 用练习3中所用观测法判断下列级数的敛散性 (1) clear;clc; epsl=0.000001; N=50000;p=1000; syms n Un=1/(n^2+n^3); s1=symsum(Un,1,N); s2=symsum(Un,1,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if sa<epsl disp('收敛') else disp('发散') end 级数1/(n^3 + n^2)收敛 clear;close syms n s=[]; for k=1:100 s(k)=symsum(1/(n^3 + n^2),1,k); end plot(s,'.') (2) clear;clc; epsl=0.000001; N=50000;p=1000; syms n Un=1/(n*2^n); s1=symsum(Un,1,N); s2=symsum(Un,1,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if sa<epsl disp('收敛') else disp('发散') end 级数1/(2^n*n)收敛 clear;close syms n s=[]; for k=1:100 s(k)=symsum(1/(2^n*n),1,k); end plot(s,'.') (3) clear;clc; epsl=0.00000000000001; N=50000;p=100; syms n Un=1/sin(n); s1=symsum(Un,1,N); s2=symsum(Un,1,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if abs(sa)<epsl disp('收敛') else disp('发散') end 级数1/sin(n)发散 clear;close syms n s=[]; for k=1:100 s(k)=symsum(1/sin(n),1,k); end plot(s,'.') 发散 (4) clear;clc; epsl=0.0000001; N=50000;p=1000; syms n Un=log(n)/(n^3); s1=symsum(Un,1,N); s2=symsum(Un,1,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if sa<epsl disp('收敛') else disp('发散') end 级数log(n)/n^3收敛 clear;close syms n s=[]; for k=1:100 s(k)=symsum(log(n)/n^3,1,k); end plot(s,'.') (5) clear;close syms n s=[];he=0; for k=1:100 he=he+factorial(k)/k^k; s(k)=he; end plot(s,'.') (6) clear;clc; epsl=0.0000001; N=50000;p=1000; syms n Un=1/log(n)^n; s1=symsum(Un,3,N); s2=symsum(Un,3,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if sa<epsl disp('收敛') else disp('发散') end 级数1/log(n)^n收敛 clear;close syms n s=[]; for k=3:100 s(k)=symsum(1/log(n)^n,3,k); end plot(s,'.') (7) clear;clc; epsl=0.0000001; N=50000;p=100; syms n Un=1/(log(n)*n); s1=symsum(Un,3,N); s2=symsum(Un,3,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if (sa)<epsl disp('收敛') else disp('发散') end 级数1/(n*log(n))发散 clear;close syms n s=[]; for k=3:300 s(k)=symsum(1/(n*log(n)),2,k); end plot(s,'.') (8) clear;clc; epsl=0.0000001; N=50000;p=100; syms n Un=(-1)^n*n/(n^2+1); s1=symsum(Un,3,N); s2=symsum(Un,3,N+p); sa=vpa(s2-s1); sa=setstr(sa); sa=str2num(sa); fprintf('级数') disp(Un) if (sa)<epsl disp('收敛') else disp('发散') end 级数((-1)^n*n)/(n^2 + 1)收敛 clear;close syms n s=[]; for k=3:300 s(k)=symsum((-1)^n*n/(n^2+1),2,k); end plot(s,'.') 实验四 积分 【练习与思考】 1.(不定积分)用int计算下列不定积分,并用diff验证 ,,,, 解:Matlab代码为: syms x y1=x*sin(x^2); y2=1/(1+cos(x)); y3=1/(exp(x)+1); y4=asin(x); y5=sec(x)^3; f1=int(y1) f2=int(y2) f3=int(y3) f4=int(y4) f5=int(y5) dy=simplify(diff([f1;f2;f3;f4;f5])) dy = x*sin(x^2) tan(x/2)^2/2 + 1/2 1/(exp(x) + 1) asin(x) (cot(pi/4 + x/2)*(tan(pi/4 + x/2)^2/2 + 1/2))/2 + 1/(2*cos(x)) + tan(x)^2/cos(x) f1 = -cos(x^2)/2 f2 = tan(x/2) f3 = x - log(exp(x) + 1) f4 = x*asin(x) + (1 - x^2)^(1/2) f5 = log(tan(pi/4 + x/2))/2 + tan(x)/(2*cos(x)) 2.(定积分)用trapz,quad,int计算下列定积分 ,,, 解:Matlab代码为 clear; x=(0+eps):0.05:1; y1=sin(x)./x; f1=trapz(x,y1) f1 =0.9460 fun1=@(x)sin(x)./x; f12=quad(fun1,0+eps,1) f12 = 0.9461 f13=vpa(int('sin(x)/x',0,1),5) f13 =0.94608 3.(椭圆的周长) 用定积分的方法计算椭圆的周长 解:椭圆的参数方程为 由参数曲线的弧长公式得 Matlab代码为 s=vpa(int('sqrt(5*sin(t)^2+4)','t',0,2*pi),5) s = 15.865 4.(二重积分)计算数值积分 解:fxy=@(x,y)1+x+y;ylow=@(x)1-sqrt(1-x.^2);yup=@(x)1+sqrt(1-x.^2); s=quad2d(fxy,-1,1,ylow,yup) s = 6.2832 或符号积分法: syms x y xi=int(1+x+y,y,1-sqrt(1-x^2),1+sqrt(1-x^2)); s=int(xi,x,-1,1) s = 2*pi 5.(假奇异积分)用trapz,quad8计算积分,会出现什么问题?分析原因,并求出正确的解。 解:Matlab代码为 clear x=-1:0.05:1; y=x.^(1/3).*cos(x); s1=trapz(x,y) fun5=@(x)x.^(1/3).*cos(x); s2=quad(fun5,-1,1) int('x^(1/3)*cos(x)','x',-1,1) s1 = 0.9036 + 0.5217i s2 = 0.9114 + 0.5262i Warning: Explicit integral could not be found. ans = int(x^(1/3)*cos(x), x = -1..1) ,原函数不存在,不能用int函数运算。 用梯形法和辛普森法计算数值积分时,由于对负数的开三次方运算结果为复数,所以导致结果错误且为复数; 显然被积函数为奇函数,在对称区间上的积分等于0,此时可以这样处理: (1)重新定义被积函数 %fun5.m function y=fun5(x) [m,n]=size(x); for k=1:m for l=1:n y(k,l)=nthroot(x(k,l),3)*cos(x(k,l)); end end end 用辛普森法: s=quad('fun5',-1,1) s = 0 用梯形法 clear; x=-1:0.01:1; y=fun5(x); s=trapz(x,y) s = -1.3878e-017 6.(假收敛现象)考虑积分, (1)用解析法求; clear; syms x k; Ik=int(abs(sin(x)),0,k*pi) Warning: Explicit integral could not be found. Ik = int(abs(sin(x)), x = 0..pi*k) (2)分别用trapz,quad和quad8求和,发现什么问题? clear; for k=4:2:8; x=0:pi/1000:k*pi; y=abs(sin(x)); trapz(x,y) end ans = 8.0000 ans = 12.0000 ans = 16.0000 for k=4:2:8 fun6=@(x)abs(sin(x)); quad(fun6,0,k*pi) end ans = 8.0000 ans = 12.0000 ans = 16.0000 7.(Simpson积分法)编制一个定步长Simpson法数值积分程序.计算公式为 其中为偶数, 解:Matlab代码为 %fun7.m function y=fun7(f_name,a,b,n) %f_name为被积函数 %[a,b]为积分区间 %n为偶数,用来确定步长h=(b-a)/n if mod(n,2)~=0 disp('n必须为偶数') return; end if nargin<4 n=100; end if nargin<3 disp('请输入积分区间') end if nargin==0 disp('error') end h=(b-a)/n; x=a:h:b; s=0; for k=1:n+1 if k==1||k==(n+1) xishu=1; elseif mod(k,2)==0 xishu=4; else xishu=2; end s=s+feval(f_name,x(k))*xishu; end y=s*h/3; end 8.(广义积分)计算广义积分 ,, 并验证公式 . 解:Matlab代码为 clear; syms x s1=vpa(int(exp(-x^2)/(1+x^4),-inf,inf),5) s2=quad(@(x)tan(x)./sqrt(x),0,1) s3=quad(@(x)sin(x)./sqrt(1-x.^2),0,1) s4=vpa(int(exp(-x^2/2)/sqrt(2*pi),-inf,inf),5) s5=int(sin(x)./x,0+eps,inf) s1 = 1.4348 s2 = 0.7968 s3 = 0.8933 s4 = 1.0 s5 = pi/2 - sinint(1/4503599627370496) 实验五 二元函数的图形 【练习与思考】 1. 画出空间曲线在范围内的图形,并画出相应的等高线。 clear; x=-30:0.5:30;y=-30:0.5:30; [X,Y]=meshgrid(x,y); Z=10*sin(sqrt(X.^2+Y.^2))./sqrt(1+X.^2+Y.^2); mesh(X,Y,Z) 2. 根据给定的参数方程,绘制下列曲面的图形。 a) 椭球面 clear; u=0:pi/50:2*pi; v=0:pi/50:pi; [U,V]=meshgrid(u,v); x=3*cos(U).*sin(V); y=2*cos(U).*cos(V); z=sin(U); mesh(x,y,z) b) 椭圆抛物面 clear; u=0:pi/50:pi/4; v=0:pi/50:2*pi; [U,V]=meshgrid(u,v); x=3*U.*sin(V); y=2*U.*cos(V); z=4*U.^2; mesh(x,y,z) axis equal c) 单叶双曲面 cle
展开阅读全文

开通  VIP会员、SVIP会员  优惠大
下载10份以上建议开通VIP会员
下载20份以上建议开通SVIP会员


开通VIP      成为共赢上传

当前位置:首页 > 教育专区 > 其他

移动网页_全站_页脚广告1

关于我们      便捷服务       自信AI       AI导航        抽奖活动

©2010-2025 宁波自信网络信息技术有限公司  版权所有

客服电话:4009-655-100  投诉/维权电话:18658249818

gongan.png浙公网安备33021202000488号   

icp.png浙ICP备2021020529号-1  |  浙B2-20240490  

关注我们 :微信公众号    抖音    微博    LOFTER 

客服