当前位置：文档库 › MATLAB实验参考答案

MATLAB实验参考答案

MATLAB程序设计与应用（第二版）

实验参考答案

%实验一MATLAB运算基础

%第一题

%（1）

z1=2*sin(85*pi/180)/(1+exp(2)) %（2）

x=[2,1+2i;-0.45,5];

z2=0.5*log(x+sqrt(1+x.^2))

%(3)

a=-3.0:0.1:3.0;

z3=(exp(0.3*a)-exp(-0.3*a))/2.* sin(a+0.3)+log((0.3+a)/2)

%(4)

t=0:0.5:2.5;

z4=t.^2.*(t>=0&t<1)+(t.^2-1).*( t>=1&t<2)+(t.^2-2*t+1).*(t>=2&t <3)

%第二题

A=[12 34 -4;34 7 87;3 65 7];

B=[1 3 -1;2 0 3;3 -2 7];

A+6*B

A-B+eye(size(A))

A*B

A.*B

A^3

A.^3

A/B

B\A

[A,B]

[A([1,3],:);B^2]

%第三题

A=[1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24 25]

B=[3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11]

C=A*B

F=size(C)

D=C(F(1)-2:F(1),F(2)-1:F(2)) whos %第四题

%(1):

A=100:999;

B=rem(A,21);

C=length(find(B==0))

%(2):

A='lsdhKSDLKklsdkl';

k=find(A>='A'&A<='Z');

A(k)=[]

%实验二MATLAB矩阵分析与处理

%第一题

E=eye(3);

R=rand(3,2);

O=zeros(2,3);

S=diag([2,3]);

A=[E,R;O,S];

A^2

B=[E,(R+R*S);O,S^2]

%第二题

H=hilb(5)

P=pascal(5)

Hh=det(H)

Hp=det(P)

Th=cond(H)

Tp=cond(P)

%第三题:

A=fix(10*rand(5))

H=det(A)

Trace=trace(A)

Rank=rank(A)

Norm=norm(A)

%第四题:

A=[-29,6,18;20,5,12;-8,8,5] [V,D]=eig(A)

%数学意义略

%第五题方法一:

%(1):

A=[1/2,1/3,1/4;1/3,1/4,1/5;1/4, 1/5,1/6];

b=[0.95,0.67,0.52]';

x=inv(A)*b

%(2):

B=[0.95,0.67,0.53]';

x=inv(A)*B

%(3):

cond(A)

%第五题方法二:

A=hilb(4)

A(:,1)=[]

A(4,:)=[]

B=[0.95,0.67,0.52]';

X=inv(A)*B

B1=[0.95,0.67,0.53]';

X1=inv(A)*B1

N=cond(B)

N1=cond(B1)

Na=cond(A) %矩阵A为病态矩阵

%第六题

A=[1,4,9;16,25,36;49,64,81]

B=sqrtm(A)

C=sqrt(A) %sqrtm函数是以矩阵为单位进行计算，sqrt函数是以矩阵中的元素进行计算

%实验三选择程序结构设计

%第一题程序一

x=[-5.0,-3.0,1.0,2.0,2.5,3.0,5. 0];

y=[];

for x0=x

if x0<0&x0~=-3

y=[y,x0*x0+x0-6];

elseif x0>=0&x0<5&x0~=2&x0~=3 y=[y,x0*x0-5*x0+6];

else

y=[y,x0*x0-x0-1];

end

x=[-5.0,-3.0,1.0,2.0,2.5,3.0,5.0];

y=[]; %建立存放所有y值的矩阵

for x0=x

if x0<0&x0~=-3

y=[y,x0*x0+x0-6];

elseif

x0>=0&x0<5&x0~=2&x0~=3

y=[y,x0*x0-5*x0+6];

else

y=[y,x0*x0-x0-1];

end

x %输出所有x

y %输出所有y

%第一题程序二

x=[-5,-3,1,2,2.5,3,5];

y=[];

for a=1:7

if x(a)<0&x(a)~=-3

y=[y,(x(a))^2+x(a)-6];

elseif

x(a)>=0&x(a)<5&x(a)~=2&x(a)~=3 y=[y,(x(a))^2-5*x(a)+6];

else

y=[y,x(a)*x(a)-x(a)-1];

end

%第二题程序一:

score=input('请输入一个分数:'); if score<0&score>100

disp('输入的分数不合理');

else

if score>=90&score<=100

disp('等级为:A');

elseif score>=80&score<=89

disp('等级为:B');

elseif score>=70&score<=79

disp('等级为:C');

elseif score>=60&score<=69

disp('等级为:D');

else

disp('等级为:E');

end

x=input('请输入一个百分制成绩：'); if x>100|x<0

disp('您输入的成绩不是百分制成绩，请重新输入。');

else

if x<=100&x>=90

disp('A');

elseif x<=89&x>=80

disp('B');

elseif x<=79&x>=70

disp('C');

elseif x<=69&x>60

disp('D');

else

disp('E');

end

%第二题程序二：

score=input('请输入一个分数:'); if score<0&score>100

disp('输入的分数不合理');

else

switch fix(score/10)

case {9,10}

disp('A');

case {8}

disp('B');

case {7}

disp('C');

case {6}

disp('D');

otherwise

disp('E');

end

end s=input('请输入一个成绩（0分到100分之间）：'); %s用于存放成绩

while

1 %判断输入成绩的合理性

if s<0|s>100

disp('输入的成绩需在0到100之间，请重新输入：')

s=input('请输入一个成绩（0分到100分之间）：');

else

break;

end

switch

fix(s/10) %对成绩做出等级判断

case {9,10}

disp('A')

case 8

disp('B')

case 7

disp('C')

case 6

disp('D')

otherwise

disp('E')

end

%第三题

a=input('工号:');

b=input('工时:');

if b>120

c=120*84+(b-120)*84*(1+0.15); elseif b<60

c=84*b-700;

else

c=84*b;

end

disp(['工号为',num2str(a),'的员工的工资为',num2str(c)]);

n=input('请输入员工工号：');

h=input('该员工工作时数是：');

if h>120

x=(h-120)*84*(1+0.15)+120*84; elseif h<60

x=h*84-700;

else

x=h*84;

end

disp([num2str(n),'号员工','的应发工资为',num2str(x)]);

%第四题（还可以用switch语句实现）

a=input('a=');

b=input('b=');

c=input('输入一个四则运算符号：','s');

if c=='+'

d=a+b;

elseif c=='-'

d=a-b;

elseif c=='*'

d=a*b;

elseif c=='/'

d=a/b;

end

disp(['结果为:',num2str(d)]);

a=fix(10+(99-10)*rand(1,2)) %产生两个随机整数

x=a(1);

y=a(2);

t=input('请输入运算符号：','s'); if t=='+'

z=x+y;

elseif t=='-'

z=x-y;

elseif t=='*'

z=x*y;

elseif t=='/'

z=x/y;

end disp([num2str(x),t,num2str(y),' =',num2str(z)]) %输出运算结果

%第五题

a=rand(5,6);

n=input('输出矩阵的第几行:');

if n>5

disp(['超出了矩阵的行数,输出矩阵的最后一行为:',num2str(a(5,:))]);

else

disp(['矩阵的第',num2str(n),'行为:',num2str(a(n,:))]);

end

a=rand(5,6) %产生5x6的随机矩阵

n=input('请输入您要输出矩阵的第几行：');

if n>5

disp('超出了矩阵的行数，矩阵的最后一行为：')

a(5,:)

else

disp(['矩阵的第',num2str(n),'行为：'])

a(n,:)

end

%实验四循环结构程序设计

%第一题程序一

s=0;

n=input('n=?');

for i=1:n

s=s+1/i/i;

end

PI=sqrt(6*s)

%第一题程序二

n=input('n=?');

a=1:n;

b=1./a.^2;

PI=sqrt(6*sum(b))

%第二题

y=0;

n=1;

while(y<3)

y=y+1/(2*n-1);

n=n+1;

end

y=y-1/(2*(n-1)-1)

n=n-2

%第三题

a=input('a=?');

b=input('b=?');

Xn=1;

Xn1=a/(b+Xn);

n=0;

while abs(Xn1-Xn)>1e-5

Xn=Xn1;

Xn1=a/(b+Xn);

n=n+1;

if n==500

break;

end

Xn1

r1=(-b+sqrt(b*b+4*a))/2

r2=(-b-sqrt(b*b+4*a))/2

%第四题

for i=1:100

if i==1

f(i)=1;

elseif i==2

f(i)=0;

elseif i==3

f(i)=1;

else

f(i)=f(i-1)-2*f(i-2)+f(i-3) ;

end

max(f)

min(f)

sum(f)

length(find(f>0))

length(find(f==0))

length(find(f<0))

%第五题：

s=0;n=0;

for i=2:49

b=i*(i+1)-1;

m=fix(sqrt(b));

for j=2:m

if rem(b,j)==0

break

end

if j==m

n=n+1;

s=s+b;

end

%实验五函数文件

%第一题

function y=mat1(x) %建立函数文件mat1.m

y=[exp(x),log(x),sin(x),cos(x)] ;

%在命令窗口调用上述函数文件：

y=mat1(1+i)

%第二题程序一

function

[a,b,N,M]=shiyanwu2(m,n,t)

A=[m*cos(t*pi/180),-m,-sin(t*pi /180),0;m*sin(t*pi/180),0,cos(t *pi/180),0;0,n,-sin(t*pi/180),0 ;0,0,-cos(t*pi/180),1];

B=[0,9.8*m,0,9.8*n];

C=inv(A)*B';

a=C(1);

b=C(2);

N=C(3);

M=C(4);

%在命令窗口调用该函数文件:

m1=input('m1=');

m2=input('m2=');

theta=input('theta=');

[a1,a2,N1,N2]=shiyanwu2(m1,m2,t heta)

%第二题程序二

function X=mat2(m1,m2,t)

g=9.8;

A=[m1*cos(t*pi/180),-m1,-sin(t* pi/180),0;m1*sin(t*pi/180),0,co s(t*pi/180),0;0,m2,-sin(t*pi/18 0),0;0,0,-cos(t*pi/180),1];

B=[0;m1*g;0;m2*g];

X=inv(A)*B;

%在命令窗口调用该函数文件：

X=mat2(1,1,60)

%第三题

function flag=mat3(x)

flag=1;

for i=2:sqrt(x)

if rem(x,i)==0

flag=0;

break;

end

%在命令窗口调用该函数文件：

for i=10:99

j=10*rem(i,10)+fix(i/10);

if mat3(i)&mat3(j)

disp(i)

end

%第四题

function y=fx(x)

y=1./((x-2).^2+0.1)+1./((x-3).^ 4+0.01);

%在命令窗口调用该函数文件：

y=fx(2)

a=[1,2;3,4];

y=fx(a) %第五题

%(1)

function f1=mat5(n)

f1=n+10*log(n*n+5);

%在命令窗口中调用该函数文件：

y=mat5(40)/(mat5(30)+mat5(20)) %(2)方法一

function f2=mat6(n)

f2=0;

for i=1:n

f2=f2+i*(i+1);

end

%在命令窗口中调用该函数文件如：

y=mat6(40)/(mat6(30)+mat6(20)) %(2)方法二

function f2=mat7(n)

i=1:n;

m=i.*(i+1);

f2=sum(m);

end

%在命令窗口中调用该函数文件如：

y=mat7(40)/(mat7(30)+mat7(20)) %实验六高层绘图操作

%第一题：

x=linspace(0,2*pi,101);

y=(0.5+3*sin(x)./(1+x.^2)).*cos (x);

plot(x,y)

%第二题：

%（1）

x=linspace(-2*pi,2*pi,100);

y1=x.^2;

y2=cos(2*x);

y3=y1.*y2;

plot(x,y1,'b-',x,y2,'r:',x,y3,' y--');

text(4,16,'\leftarrow y1=x^2'); text(6*pi/4,-1,'\downarrow

y2=cos(2*x)');

text(-1.5*pi,-2.25*pi*pi,'\upar row y3=y1*y2');

%（2）

x=linspace(-2*pi,2*pi,100);

y1=x.^2;

y2=cos(2*x);

y3=y1.*y2;

subplot(1,3,1);%分区

plot(x,y1);

title('y1=x^2');%设置标题

subplot(1,3,2);

plot(x,y2);

title('y2=cos(2*x)');

subplot(1,3,3);

plot(x,y3);

title('y3=x^2*cos(2*x)');

%（3）

x=linspace(-2*pi,2*pi,20);

y1=x.^2;

subplot(2,2,1);%分区

bar(x,y1);

title('y1=x^2的条形图');%设置标题subplot(2,2,2);

stairs(x,y1);

title('y1=x^2的阶梯图');

subplot(2,2,3);

stem(x,y1);

title('y1=x^2的杆图');

subplot(2,2,4);

fill(x,y1,'r');%如果少了'r'则会出错

title('y1=x^2的填充图');

%其他的函数照样做。

%第三题

x=-5:0.01:5;

y=[];%起始设y为空向量

for x0=x

if x0<=0 %不能写成x0=<0

y=[y,(x0+sqrt(pi))/exp(2)]; %将x对应的函数值放到y中

else

y=[y,0.5*log(x0+sqrt(1+x0^2 ))];

end

plot(x,y)

%第四题：

a=input('a='); b=input('b=');

n=input('n=');

t=-2*pi:0.01:2*pi;

r=a*sin(b+n*t);

polar(t,r)

%第五题

x=linspace(-5,5,21);

y=linspace(0,10,31);

[x,y]=meshgrid(x,y);%在[-5,5]*[0,10]的范围内生成网格坐标

z=cos(x).*cos(y).*exp(-sqrt(x.^ 2+y.^2)/4);

subplot(2,1,1);

surf(x,y,z);

subplot(2,1,2);

contour3(x,y,z,50);%其中50为高度的等级数，越大越密

%第六题

ezsurf('cos(s)*cos(t)','cos(s)* sin(t)','sin(s)',[0,0.5*pi,0,1. 5*pi]); %利用ezsurf隐函数

shading interp %进行插值着色处理%实验七低层绘图操作

%第一题

h=figure('MenuBar','figure','co lor','r','WindowButtonDownFcn', 'disp(''Left Button Pressed'')') %第二题

x=-2:0.01:2;

y=x.^2.*exp(2*x);

h=line(x,y);

set(h,'color','r','linestyle',' :','linewidth',2)

text(1,exp(2),'y=x^2*exp(2*x)') %第三题

t=0:0.00001:0.001;

[t,x]=meshgrid(t);

v=10*exp(-0.01*x).*sin(2000*pi* t-0.2*x+pi);

axes('view',[-37.5,30]);

h=surface(t,x,v);

title('v=10*exp(-0.01*x).*sin(2 000*pi*t-0.2*x+pi)');

xlabel(Ct'),ylabel('x'),zlabel( 'v')

%第四题

x=0:0.01:2*pi;

y1=sin(x);

y2=cos(x);

y3=tan(x);

y4=cot(x);

subplot(2,2,1);

plot(x,y1);

subplot(2,2,2);

plot(x,y2);

subplot(2,2,3);

plot(x,y3);

subplot(2,2,4);

plot(x,y4);

%第五题

cylinder(5);

light('Position',[0,1,1]); material shiny

%实验八数据处理与多项式运算

%第一题

%(1)

A=rand(1,30000);

b=mean(A)

std(A,0,2)

%(2)

max(A)

min(A)

%(3)

n=0;

for i=1:30000

if A(i)>0.5

n=n+1;

end

p=n/30000

%第二题

%(1)

A=45+51*rand(100,5);

[Y,U]=max(A) [a,b]=min(A)

%(2)

m=mean(A)

s=std(A)

%(3)

sum(A,2)

[Y,U]=max(ans)

[a,b]=min(ans)

%(4)

[zcj,xsxh]=sort(ans)

%第三题

h=6:2:18;

x=6.5:2:17.5;

t1=[18,20,22,25,30,28,24];

t2=[15,19,24,28,34,32,30];

T1=spline(h,t1,x)

T2=spline(h,t2,x)

%第四题

x=1:0.1:101;

y1=log10(x);

p=polyfit(x,y1,5)

y2=polyval(p,x);

plot(x,y1,':',x,y2,'-')

%第五题

%(1)

p1=[1,2,4,0,5];

p2=[1,2];

p3=[1,2,3];

p=p1+[0,conv(p2,p3)] %为使两向量大小相同，所以补0

%(2)

A=roots(p)

%(3)

A=[-1,1.2,-1.4;0.75,2,3.5;0,5,2 .5];

polyval(p,A)

%(4)

polyvalm(p,A)

实验十

程序：

x=sym('6');

y=sym('5');

z=(x+1)/(sqrt(3+x)-sqrt(y))

1、分解因式

（1）

程序：

syms x y;

A=x^4-y^4;

factor(A)

（2）

程序：

factor(sym('5135'))

3、化简表达式

（1）

程序：

syms beta1 beta2

y=sin(beta1)*cos(beta2)-cos(beta 1)*sin(beta2)

simple(y)

（2）

程序：

syms x

y=(4*x^2+8*x+3)/(2*x+1)

simple(y)

5、用符号方法求下列极限或导数

（1）

程序：

syms x

f=(x*(exp(sin(x))+1)-2*(exp(tan(x)) -1))/(sin(x))

limit(f)

（2）

程序：

syms x

y=(sqrt(pi)-sqrt(acos(x)))/(sqrt(x+1 ));

limit(f,x,-1,'right')

（3）

程序：

syms x

y=(1-cos(2*x))/x;

y1=diff(y)

y2=diff(y,x,2)

6、用符号方法求下列积分

（1）

程序：

syms x f=1/(1+x^4+x^8)

int(f)

（2）

程序：

syms x

f=1/(((asin(x))^2)*sqrt(1-x^2))

int(f)

（3）

程序：

syms x

f=(x^2+1)/(x^4+1)

int(f,x,0,inf)

（4）

程序：

syms x

f=exp(x)*(1+exp(x))^2

y=int(f,x,0,log(2))

double(y)

实验十一级数与方程符号求解

1. 级数符号求和。

(1) 计算。

(2) 求级数的和函数，并求之和。

解：

M文件：

clear all;clc;

n=sym('n');x=sym('x');

S1=symsum(1/(2*n-1),n,1,10)

S2=symsum(n^2*x^(n-1),n,1,inf)

S3=symsum(n^2/5^n,n,1,inf) %vpa(S3)可以转化成小数

运行结果：

S1 =

31037876/14549535

S2 =

piecewise([abs(x) < 1, -(x^2 + x)/(x*(x - 1)^3)])

S3 =

15/32

2. 将lnx在x=1处按5次多项式展开为泰勒级数。

解：

M文件：

clear all;clc;

x=sym('x');

f=log(x);

taylor(f,x,6,1)

运行结果：

ans =

x - (x - 1)^2/2 + (x - 1)^3/3 - (x - 1)^4/4 + (x - 1)^5/5 - 1

3. 求下列方程的符号解。

解：

M文件：

clear all;clc;

x1=solve('log(x+1)-5/(1+sin(x))=2')

x2=solve('x^2+9*sqrt(x+1)-1')

x3=solve('3*x*exp(x)+5*sin(x)-78.5')

[x4

y4]=solve('sqrt(x^2+y^2)-100','3*x+5*y-8')

运行结果：

x1 =

521.67926389905839979437366649258

x2 =

-1

(3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) - (1/2*6465^(1/2) + 2171/54)^(1/3)))/2 - (6465^(1/2)/2 + 2171/54)^(1/3)/2 - 2/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) + 1/3

1/3 - (6465^(1/2)/2 + 2171/54)^(1/3)/2 - (3^(1/2)*i*(4/(9*(6465^(1/2)/2 + 2171/54)^(1/3)) - (1/2*6465^(1/2) + 2171/54)^(1/3)))/2 - 2/(9*(6465^(1/2)/2 + 2171/54)^(1/3))

x3 =

2.3599419584772910151699327715486

x4 =

12/17 - (10*21246^(1/2))/17

(10*21246^(1/2))/17 + 12/17

y4 =

(6*21246^(1/2))/17 + 20/17

20/17 - (6*21246^(1/2))/17

4. 求微分方程初值问题的符号解，并与数值解进行比较。

解：

M文件：clear all;clc;

dsolve('D2y+4*Dy+29*y','y(0)=0','Dy(0)=15',' x')

运行结果：

ans =

(3*sin(5*x))/exp(2*x)

5. 求微分方程组的通解。

解：

M文件：

clear all;clc;

[x y z]=dsolve('Dx=2*x-3*y+3*z',...

'Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t') 运行结果：

x =

C1/exp(t) + C2*exp(2*t)

y =

C1/exp(t) + C2*exp(2*t) + C3/exp(2*t)

z =

C2*exp(2*t) + C3/exp(2*t)

实验九数值微积分与方程数值求解

1. 求函数在指定点的数值导数。

解：M文件：

clc;clear;

x=1;

i=1;

f=inline('det([x x^2 x^3;1 2*x 3*x^2;0 2 6*x])');

while x<=3.01

g(i)=f(x);

i=i+1;

x=x+0.01; %以0.01的步长增加，可再缩小步长提高精度

end

t=1:0.01:3.01;

dx=diff(g)/0.01; %差分法近似求导

f1=dx(1) %x=1的数值倒数

f2=dx(101) %x=2的数值倒数

f3=dx(length(g)-1) %x=3的数值倒数

运行结果：

f1 =

6.0602

f2 =

24.1202

f3 =

54.1802

2. 用数值方法求定积分。

(1) 的近似值。

(2)

解：M文件：

clc;clear;

f=inline('sqrt(cos(t.^2)+4*sin(2*t).^2+1)');

I1=quad(f,0,2*pi)

g=inline('log(1+x)./(1+x.^2)');

I2=quad(g,0,2*pi)

运行结果：

3. 分别用3种不同的数值方法解线性方程组。

解：M文件：

clc;clear;

A=[6 5 -2 5;9 -1 4 -1;3 4 2 -2;3 -9 0 2];

b=[-4 13 1 11]';

x=A\b

y=inv(A)*b

[L,U]=lu(A);

z=U\(L\b)

运行结果：

4. 求非齐次线性方程组的通解。

解：M文件

function [x,y]=line_solution(A,b)

[m,n]=size(A);

y=[ ];

if norm(b)>0 %非齐次方程组

if rank(A)==rank([A,b])

if rank(A)==n

disp('有唯一解x');

x=A\b;

else

disp('有无穷个解，特解x，基础解系y');

x=A\b;

y=null(A,'r');

end

else

disp('无解');

x=[ ];

end

else %齐次方程组

disp('有零解x');

x=zeros(n,1);

if rank(A)

disp('有无穷个解，基础解系y');

y=null(A,'r');

end

clc;clear;

format rat

A=[2 7 3 1;3 5 2 2;9 4 1 7];

b=[6 4 2]';

[x,y]=line_solution(A,b)

运行结果：

有无穷个解，特解x，基础解系y Warning: Rank deficient, rank = 2, tol = 8.6112e-015.

> In line_solution at 11

x =

-2/11

10/11

y =

1/11 -9/11

-5/11 1/11

1 0

0 1

所以原方程组的通解是：

，其中为任意常数。

5. 求代数方程的数值解。

(1) 3x+sinx-ex=0在x0=1.5附近的根。

(2) 在给定的初值x0=1，y0=1，z0=1下，求方程组的数值解。

解：M文件：

function g=f(x)

g=3*x+sin(x)-exp(x);

clc;clear;

fzero('f',1.5)

结果是：

ans =

1289/682

(2). M文件：

function F=fun(X)

x=X(1);

y=X(2);

z=X(3);

F(1)=sin(x)+y^2+log(z)-7;

F(2)=3*x+2-z^3+1;

F(3)=x+y+z-5;

X=fsolve('myfun',[1,1,1]',optimset('Display','o ff'))

运行结果：

6. 求函数在指定区间的极值。

(1) 在(0,1)内的最小值。

(2) 在[0,0]附近的最小值点和最小值。解：M文件：

function f=g(u)

x=u(1); y=u(2);

f=2*x.^3+4*x.*y^3-10*x.*y+y.^2;

clc;clear;

format long

f=inline('(x^3+cos(x)+x*log(x))/exp(x)'); [x,fmin1]=fminbnd(f,0,1)

[U,fmin2]=fminsearch('g',[0,0])

运行结果

7. 求微分方程的数值解。

解：M文件：

function xdot= sys( x,y)

xdot=[y(2);(5*y(2)-y(1))/x];

clc;clear;

x0=1.0e-9;xf=20;

[x,y]=ode45('sys',[x0,xf],[0 0]);

[x,y]

运行结果：

8. 求微分方程组的数值解，并绘制解的曲线。

解：令y1=x,y2=y,y3=z; 这样方程变为: ,自变量是t

M文件：

function xdot=sys(x,y)

xdot=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]; clc;clear;

t0=0;tf=8; [x,y]=ode23('sys',[t0,tf],[0,1,1]) plot(x,y)