MATLAB数学实验A

clc;
a=1;b=1;
ezplot(sprintf('x^2/%f-y^2/%f',a^2,b^2));
hold on;
ezplot(sprintf('x^2/%f-y^2/%f-1',a^2,b^2));
ezsurf('sin(a)*cos(b)','sin(a)*sin(b)','cos(a)',[0,pi,0,2*pi],60); hold on;
ezsurf('x^2+y^2',[-1,1,-1,1],60);
clear all;
x=-8:0.1:8;
y=-8:0.1:8;
[X,Y]=meshgrid(x,y);
Z=sin(sqrt(X.^2+Y.^2))./sqrt(X.^2+Y.^2+2); [X,Y,Z]=peaks(50);
surf(X,Y,Z)
syms x y;
y=2*x^3-6*x^2-18*x+7; solve(diff(y,x),x)
x=-1;eval(y)
x=3;eval(y)
syms x y;
z='x*y';
dblquad(z,1,4,-1,2)
结果
ans =
11.2500
求函数1+x -exp(2*x)+5的原函数clear all
syms x C;
f=int(1+x -exp(2*x)+5,'x')+C
syms x y;
>> x=0:0.01:1;
>> y=sin(sin(x));
>> trapz(x,y)
x=0:0.05:1;
y=[1.97687 2.17002 2.34158 2.46389 2.71512 3.06045 3.27829 3.51992 3.8215 4.2435 4.55188 4.88753 5.15594 5.698 6.04606 6.42701 7.00342 7.50192 7.89178 8.49315 9.0938]
cftool
解常微分方程y’=-0.9y/(1+2x)的数值解y(0)=1
从0 到0. 1 的数值解，取步长0.02
clear all
x1=0;
x2=0.1;
h=0.02;
y(1,1)=1;
for i=1:(x2-x1)/h
y(i+1,1)=y(i,1)-0.9*h*y(i,1)/(1+2*0.02*(i-1)); end
clear
>> s=0;
>> n=1000000;
>> digits(22);
>> for k=1:n
s=s+(-1)^(k-1)/(2*k-1);
end
>> s=vpa(4*s,20)
clear all
P=((1/16)^0)*((4/(8*0+1))-(2/(8*0+4))-(1/(8*0+5))-(1/(8*0+6))); for n=1:5
P=P+((1/16)^n)*((4/(8*n+1))-(2/(8*n+4))-(1/(8*n+5))-(1/(8*n+6))); end
vpa(P,20)
absolute=abs(P-pi)
结果
ans =
3.1415926532280877836
absolute =
3.6171e-010
判断n^2+n+41所求是否有素数
clear all
n=0:100;
V=power(n,2)+n+41
P=isprime(V')
明显可以给出素数
判断在10000之内，由n^2+n+41算得的素数有多少clear all
n=0:100;
V=power(n,2)+n+41
P=isprime(V')
j=0;
for i=1:101
if(V(1,i)<=10000)
if (P(i,1)==1)
j=j+1;
else
j=j;
end
end
end
j
结果为86
判断在10000之内，由n^2-79*n+1601算得的素数有多少clear all
n=0:100;
V=power(n,2)-79*n+1601;
P=isprime(V');
j=0;
for i=1:101
if(V(1,i)<=10000)
if (P(i,1)==1)
j=j+1;
else
j=j;
end
end
end
j
结果为96
判断在10000之内，由6*n^2+6*n+31算得的素数有多少m文件：
function[d] =sushu(n)
d=0;
for i=0:n
a=0;
b=6*i.^2+6*i+31;
for c=2:b-1
if mod(b,c)==0
a=a+1;
end
end
if(a==0)
d=d+1; end
end
主程序：
>> d=sushu(100)
d =
76
简单迭代法：
a 的m 次方根m a 的迭代格式为：])1[(1111---+-=m n n n a a a m m a 程序：
m 文件：
function[b]=diedai(a,c) b=10;
while b-a>c
b=(2*a+2\(a.^2))\3; end
主程序：
>> a=1.5;
>> c=10.^-6;
>> b=diedai(a,c)
初值为1.5
二分法
x1=1;
x2=2;
e=0.0001;
y=0;
while(abs(y-2)>e)
y1=(x1+x2)/2;
y=y2^3;
if(y>2)
x2=y1;
else
x1=y1;
end
end
y1=vpa(y1,10)
结果为
y1=
1.1259918213
(1)
n=8;
x=cosd([90 210 330]');
y=sind([90 210 330]');
tri=[1 2 3];
tnum=1;dnum=3;
for ii=2:n
tri2=zeros(3*tnum,3);
nn=1;
for jj=1:tnum
xnew=(x(tri(jj,[1 2 3]))+x(tri(jj,[2 3 1])))/2;
ynew=(y(tri(jj,[1 2 3]))+y(tri(jj,[2 3 1])))/2;
x=[x;xnew];
y=[y;ynew];
tri2(nn,:)=[tri(jj,1),dnum+1,dnum+3];
tri2(nn+1,:)=[dnum+1,tri(jj,2),dnum+2];
tri2(nn+2,:)=[dnum+3,dnum+2,tri(jj,3)];
nn=nn+3;
dnum=dnum+3;
end
tri=tri2;
tnum=3*tnum;
end
triplot(tri,x,y);
axis image;
axis off;
(2)
x=[0.5;0.5];
plot(x(1),x(2),'.')
h=[0.75 0.90 0.98 1.00];
b1=[0;1.5];
b2=[0;1.5];
b3=[0;0.4];
A1=[0.75 0.04 ;-0.04 0.75];A2=[0.20 -0.3;0.25 0.23]; A3=[-0.15 0.32;0.26 0.24];A4=[0 0 ;0 0.15];
for i=1:20000
r=rand;
if r<h(1)
x=A1*x+b1;
elseif r<h(2)
x=A2*x+b2;
elseif r<h(3)
x=A3*x+b3;
else
x=A4*x;
end
plot(x(1),x(2),'g'),hold on
end
clear all
a=0.04;
for j=1:(4/0.04)
x(1,j)=0.5;
for i=1:98
x(i+1,j)=a*j*x(i,j)*(1-x(i,j));
end
end
clear i j
for j=1:100
for i=1:50
y(i,j)=x(i+49,j);
end
end
t=0.04:0.04:4;
plot(t,y,'b.')
xlabel('a');ylabel('x');title('图1')。