3.代码：
function B=chuli(A)
[m,n]=size(A);
C=[];
X=zeros(1,n);
for k=1:m
if sum(A(k,:)~=X)
C=[C;A(k,:)];
end
end
[m,n]=size(C);
B=[];
X=zeros(m,1);
for k=1:n
if sum(C(:,k)~=X)
B=[B,C(:,k)];
end
end
例：
>> A=[1 2 3 0 4;5 6 7 0 8;0 0 0 0 0;1 3 2 0 1]
A =
1 2 3 0 4
5 6 7 0 8
0 0 0 0 0
1 3 2 0 1
>> D=chuli(A)
D =
1 2 3 4
5 6 7 8
1 3 2 1
6.代码：
function y=fun6(x)
if max(x)>3 | min(x)<-3
error('定义域为[-3,3]！');
end
y=x;
n=length(x);
for k=1:n
if x(k)<-1
y(k)=(-x(k)^2-4*x(k)-3)/2;
elseif x(k)<1
y(k)=-x(k)^2+1;
else
y(k)=(-x(k)^2+4*x(k)-3)/2;
end
end
例：
>> x=-3:0.05:3;
>> y=fun6(x);
>> plot(x,y)
7. 代码：
function [y,a]=fun7(n)
a=zeros(1,n+1);
a(1)=1;a(2)=2;
y=2;
for k=3:n+1
a(k)=a(k-1)+a(k-2);
y=y+a(k)/a(k-1);
end
例：
>> [sum,a]=fun7(15)
sum =
24.570090982808093
a =
Columns 1 through 8
1 2 3 5 8 13 21 34
Columns 9 through 16
55 89 144 233 377 610 987 1597
11.代码：
%保存文件名t11.m
t=[1 2 3 4 5 6 7 8 9 10];
y=[4.842 4.362 3.754 3.368......
3.169 3.038 3.034 3.016 3.012 3.005];
tt=1:0.05:10;
u=exp(-t);
p=polyfit(u,y,1);
uu=exp(-tt);
yy1=polyval(p,uu);
z1=polyval(p,u);
wucha1=sqrt(sum((z1-y).^2))
v=t.*u;
q=polyfit(v,y,1);
vv=tt.*uu;
yy2=polyval(q,vv);
z2=polyval(q,v);
wucha2=sqrt(sum((z2-y).^2))
figure(1);
plot(t,y,'+',tt,yy1,t,z1,'o');
figure(2);
plot(t,y,'+',tt,yy2,t,z2,'o');
例：
>> t11
wucha1 =
0.727954807862586
wucha2 =
0.037478017693910
14.代码：
%保存文件名t14.m
t=[0 0.3 0.8 1.1 1.6 2.3];
y=[0.5 0.82 1.14 1.25 1.35 1.41];
tt=0:0.01:2.3;
a=polyfit(t,y,2) %系数a，从高次到低次系数
yy1=polyval(a,tt);
z1=polyval(a,t);
wucha1=sqrt(sum((z1-y).^2))
B=[ones(size(t')) exp(-t)' ( t.*exp(-t))'];
%超定方程组系数矩阵
b=B\y' %系数b
yy2=b(1)+b(2)*exp(-tt)+b(3)*tt.*exp(-tt);
z2=b(1)+b(2)*exp(-t)+b(3)*t.*exp(-t);
wucha2=sqrt(sum((z2-y).^2))
figure(1);
plot(t,y,'+',tt,yy1,t,z1,'o');
figure(2);
plot(t,y,'+',tt,yy2,t,z2,'o');
结果：
>> t14
a =
-0.234631588556964 0.913352704820586 0.532598625427850
wucha1 =
0.071996179795285
b =
1.413231229055381
-0.913632775307342
0.380809287598076
wucha2 =
0.002645208919318
15.代码：
%Van Der Pol微分方程
function dy=vdpol(t,x)
global u; % 全局变量声明，可以在函数外赋值
dy=[x(2);u*(1-x(1)^2)*x(2)-x(1)];
%操作脚本，保存文件名t15.m
global u; % 全局变量声明
u=1; % 全局变量赋值
[t,y]=ode45(@vdpol,[0,20],[2;0]);
figure(1);
plot(t,y(:,1),t,y(:,2),'-.');
legend('解图像','导函数');
title('Van Der Pol方程解u=1');
grid on % 显示坐标网格
u=5;
[t,y]=ode45(@vdpol,[0 20],[2;0]);
figure(2);
plot(t,y(:,1));
title('Van Der Pol方程解（ode45）u=5');
grid on % 显示坐标网格
u=1000;
[t,y]=ode15s(@vdpol,[0,3000],[2;0]);
figure(3);
plot(t,y(:,1));
title('Van Der Pol方程解（ode15s）u=1000');
grid on % 显示坐标网格
结果：
>>t15
注：从u=1、5到1000的过程大家可以看到解函数的图形在急剧地变得陡峭。此时导函数y’的值变得异常的大（相对于函数值y）。在解微分方程组时（二阶微分方程是转化为微分方程组的），两个分量的函数值数量级相差太大会导致误差急剧膨胀（就像大象看蚂蚁），此类方程（组）称为刚性微分方程（组）。此时经典的ode45（即，经典的4-5阶Runge-Kutta算法）失去效用，应转用针对刚性方程的ode15s或ode23等。
16.代码：
函数文件：
function y=fun16(x)
y=exp(-x)-1.5*exp(2*cos(2*pi*x));
function y=fun16fu(x)
y=-exp(-x)+1.5*exp(2*cos(2*pi*x));
操作脚本：
%保存文件名t16.m
x=-1:.01:1;
y=fun16(x);
%plot(x,y) %观察函数图像
rt=[fzero(@fun16,-0.8) fzero(@fun16,-0.2).....
fzero(@fun16,0.2) fzero(@fun16,0.6)]
rty=fun16(rt) %验证是否为0
jx=fminbnd(@fun16,-0.2,0.2)
jd=[fminbnd(@fun16fu,-0.8,-0.4) .......
fminbnd(@fun16fu,0.2,0.8)]
plot(x,y,[-1,1],[0 0],rt,fun16(rt),'o',.....
jx,fun16(jx),'*k',jd,fun16(jd),'dr');
legend('函数','x=0','根','极小点','极大点','Location','SouthWest')
结果：
>> t16
rt =
-0.779979924312218 -0.261468786317026 0.308059714503363 0.660532508241816
rty =
1.0e-015 *
-0.888178419700125 0 0 0.222044604925031
jx =
0.001129024320210
jd =
-0.587828160360702 0.462479013607980
评论