文档库

最新最全的文档下载
当前位置:文档库 > ch6第六章习题

ch6第六章习题

%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为y(0)=1,[0 1 2 3]中0不可缺

%Excercise 1(3)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t)
%运行下列指令
clear;close;
fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)];
[t,y]=ode45(fun,[0 5],[0;1]);
plot(t,y(:,1))

%Excercise 1(5)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1)
%运行下列指令,注意参数mu的处理
clear;close;
fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)];
[t,y]=ode45(fun,[0 20],[2;0],[],1);
plot(y(:,1),y(:,2));hold on;
[t,y]=ode45(fun,[0 20],[2;0],[],2);
plot(y(:,1),y(:,2),'r');hold off;

%Excercise 2
roots([1 10 54 132 137 50])
%通解A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t)

%Excercise 3
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)
%所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)
%所用节点很少


%Excercise 4
clear;
dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3; x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3]','t','x');
[t,x]=ode45(dfun,[0 24],[1.2; 0; 0; -1.04935371]);
plot(x(:,1),x(:,3));

%Excercise 5
%方程y'=2x+y^2,y(0)=0
clear;close;
fun=inline('2*x+y^2','x','y');
[x,y]=ode45(fun,[0 1.57],0); %x的上界再增加,解会"爆炸"
plot(x,y)


%Excercise 6
clear;close;
fun=@(t,x,a,b)a*x+b;
[t,x]=ode45(fun,[0 10],0.1,[],1,1);
subplot(2,4,1);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,1);
subplot(2,4,2);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],1,-1);
subplot(2,4,3);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,-1);
subplot(2,4,4);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,1);
subplot(2,4,5);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,1);
subplot(2,4,6);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,-1);
subplot(2,4,7);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,-1);
subplot(2,4,8);plot(t,x)

%Excercise 7
%微分方程 T'=k(c-T),T(0)=20
dsolve('DT=k*(c-T)','T(0)=20','t')
%得c+exp(-k*t)*(-c+20)
%利用T(10)=25.2, T(20)=28.32拟合(或者解非线性方程)
fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t')
lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])
%解得户外温度c=33,比例系数k=0.05.

%Excercise 8
%微分方程 x'/x=0.5*(1-x),x(0)=0.05
fun=inline('0.5*(1-x)*x','t','x');
[t,x]=ode45(fun,[0 10],0.05);
plot(t,x)
id=min(find(x>0.8));
t(id)

%Excercise 9
%微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t)
%答案(1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6

ch6第六章习题

(共1页)