二、已知速度曲线v(t) 上的四个数据点下表所示：1：分析：该题提示用三次样条差值求，即用spline，所求位移是对函数的积分，所以可以用分割求和法做，把时间分成301个小段，每段长度为0.0001，用spline （t，v，x）进行拟合，每个x对应一个y值，得到的数集用y表示。
用for循环计算出分割后，拟合后的y值的和，然后乘以时间间隔0.0001得出位移。
求解加速度时取时间点0.18和之前的那一时间点，用速度的差除以时间段0.0001，得出t=0.18时的加速度a。
2：代码如下:t=[0.15:0.01:0.18];v=[3.5 1.5 2.5 2.8];t0=(0.15:0.0001:0.18);v0=spline(t,v,t0);sum=0;for j=1:1:301sum=v0(j)+sum;ends=sum*0.0001a=(v0(301)-v0(300))/0.0001plot(t0,v0)hold on;plot(t0(1),v0(1),'r*')hold on;plot(t0(101),v0(101),'r*')hold on;plot(t0(201),v0(201),'r*')hold on;plot(t0(301),v0(301),'r*')3：运行结果截图：显示结果为：s =0.0689a =-126.13954：做题中遇到的困难和收获：开始对三次样条插值函数的使用不熟悉，不明白其中的的含义，后来通过翻看matlab参考书解决了遇到的问题。
同时知道了求导求积分的一种方法，通过这次编程，我对插值，拟合的原理有所了解，学会使用三次样条插值spline（）函数并且，复习了for循环，了解了matlab的更多的功能。
在这次编程中，开始由于分段太少而导致图像曲线不圆滑，经过改进，得到了上述图像，求加速度的时候由于没看清题中是求t=0.18时的加速度，而导致加速度用变量表示，没能正确求得，经改正得到上述加速度。
三、计算图片文件 tu.bmp 给出的两个圆 A，B 的圆心，和两个圆的两条外公切线和两条内公切线的切点的坐标。
1：分析：(1) 先求圆心，调用并显示两个圆的图像，把圆用矩阵分割，转换成double数据类型，用mean（）函数求圆心，代码如下：data=imread('E:\tu.bmp')data2=double(data(1:200,1:200))data3=double(data(150:350,250:500))imshow('E:\tu.bmp')[x1,y1]=find(data2==4)%第一个圆x1center=mean(x1)%mean是求均值y1center=mean(y1)hold onplot(y1center,x1center,'bo')hold on[x2,y2]=find(data3==4)x2center=mean(x2)+150y2center=mean(y2)+250plot(y2center,x2center,'bo')(2)然后求圆的半径，（两条相互垂直的直径相加然后求均值）代码如下：%求圆的半径；r1=(max(x1)-min(x1)+max(y1)-min(y1))/4r2=(max(x2)-min(x2)+max(y2)-min(y2))/4(3)利用几何关系导出切线与圆心坐标（点到直线的距离）及半径之间的关系，求得切线方程，用solve()解方程组及eval()语句编写，（注意x和y的顺序）。
代码如下：%利用圆心到切线的距离为半径，解方程组FC,得到四条直线；FC=solve('(a*y1center-x1center+b)^2-(1+a^2)*r1^2','(a*y2center-x2ce nter+b)^2-(1+a^2)*r2^2','a','b');a=eval(FC.a)%a,b均为矩阵，为两圆的四条公切线b=eval(FC.b)%执行FC对a和b赋值x=0:0.5:530;L1=a(1)*x+b(1);L2=a(1)*x+b(2);L3=a(3)*x+b(3);L4=a(4)*x+b(4);hold onplot(x,L1,'g-',x,L2,'g-',x,L3,'g-',x,L4,'g-')（4）最后，由切线及圆心坐标和半径，根据几何关系，导出切点坐标的计算公式，利用两次for语句循环，分别求出两个圆的四个切点坐标：代码如下：%求切点坐标；%第一个圆的切点QD；for i=1:4%一个圆上四个切点，用循环求解QD=solve('(x-y1center)^2+(y-x1center)^2-r1^2','a(i)*x+b(i)-y','x',' y');x1i=eval(QD.x)%切点坐标y1i=eval(QD.y)plot(x1i,y1i,'*');end%第二个圆的切点QD，原理同第一个。
for i=1:4QD=solve('(x-y2center)^2+(y-x2center)^2-r1^2','a(i)*x+b(i)-y','x',' y');x2i=eval(QD.x)y2i=eval(QD.y)plot(x2i,y2i,'*');end2：运行结果截图：A的半径：r1 =79.7500；B的半径：r2 = 80；切线方程：y=-0.0033x+166.5824；y= 2.8964x+-475.2352；y= 0.7088x+ -88.9399；y= 0.7116x+ 106.3913；A的切点：（109．82，166.22）、（184.94，60.44）、（155.68 ,21.40）、（63.32，151.45）；B的切点：（333.20，165.48）、（257.84,271.59）、（379.72,180.22）、（287.08,310.67）3：做题中遇到的困难和收获：（1）在编程过程中很多语句不会使用，例如求切线方程组的solve函数等，在求切点坐标时，由于不熟悉for语句而一个个求，造成了很多麻烦，经过修改得上述代码。
经常搞错x和y的顺序，所以程序开始运行时比较混乱，经分析捋顺最后得到结果；(2)学会了使用solve()函数及eval()函数；(3)对求圆心之类的问题有了大体的了解，掌握了更多有关matlab的知识。
四、求微分方程组的数值解，并画出解曲线1：分析：化简方程设 y(1)=x(t) dy(1)=x’(t)=dx/dty(2)=y(t) dy(2)=y’(t)=dy/dty(3)=z(t) dy(3)=z’(t)=dz/dt则方程可化为：dy(1)=-10*y(1)+10*y(2);dy(2)=28*y(1)-y(2)-y(1)*y(3);dy(3)=-(8/3)*y(3)+y(1)*y(2);2：建一个M文件，定义方程函数，注意文件名和函数名要相同。
代码如下:function dy=vdp(t,y);dy=zeros(3,1);dy(1)=-10*y(1)+10*y(2);dy(2)=28*y(1)-y(2)-y(1)*y(3);dy(3)=-(8/3)*y(3)+y(1)*y(2);然后在命令窗口中运行：t=[0,30];% t的初值和末值y0=[1,0,0];[t,y]=ode45('vdp',t,y0);plot(t,y(:,1),'b',t,y(:,2),'g',t,y(:,3),'r')3:运行结果截图：4：做题中遇到的困难和收获：在刚开始把问题分析出来写M文件时，以为题中含有三个变量和一个中间值t,就要把图画成三维立体的，为了最后一行画图的代码，差点把整个程序都放弃，经过反复修改，最后把图画在了二维坐标系中得到了如上结果。
在调试过程中，犯了很多低级错误，例如，把方程抄错了，得到的图与预期的不一样，等等，由于有老师写的例题，所以在写程序方面有一定参照，少走了很多弯路。
通过这次编程我学会了对函数的定义，用matlab 解方程及解方程的函数命令ode45（）等。
五、预测2012-2020年威海人口数量。
说明数据的可靠来源，给出模型的假设，模型及求解的代码和计算结果。
想想如何验证你模型的准确性。
1：分析：模型假设：假设威海市人口的变化规律满足Logistic人口预测模型。
模型建立：设人口数量为x；人口平均增长率为r;人口固有增长率为r（x）是与x有关的函数；最大人口数量（人口容量）为xm；dx/dt=r(x)x,x(0)=x0r(x)=r-sxx=xm时人（xm）=0，于是s=r/xm,r(x)=r(1-x/xm)dx/dt=rx(1-x/xm),x(0)=x0;最后得到：x(t)=xm/(1+(xm/x0-1)e^(-rt)).2:建一个M文件进行人口拟合，前半部分是拟合和检测，后半部分是预测data = [2189973 2213869 2238989 2251553 2258542 2270689 2301232 23257392342838 2361453 2374563 2385968 2393181 2406062 2419678 2430043 2439975 2450927 2457537 2462245 2469542 2472215 2476178 2476269 2483889 2490904 2498299 2510553 2522307 2529677];time = [1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009];plot(time,data,'--');pData = data(1:28);%前二十八个数拟合，后两个数检测t = linspace(1,30,30);dData = pData(2:end) - pData(1:end-1);%后一年与前一年的差值y = dData./pData(2:end);%增长率p = polyfit(pData(2:end),y,1);%拟合，1是拟合的次数，得到地p有两个数p1，p2分别是参数s=-p(1);r = p(2);xm = r/s;data2 = xm./(1+(xm./pData(1)-1).*exp(-r*t))%pdata(1)即x的初值hold on;plot(time,data2) %对2010-2020年数据进行拟合tm2=[1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020];pData =data(1:30);%前三十个数的拟合t = linspace(1,41,41);%对四十个数进行拟合dData = pData(2:end) - pData(1:end-1);y = dData./pData(2:end);p = polyfit(pData(2:end),y,1);s=-p(1);r = p(2);xm = r/s;data2 = xm./(1+(xm./pData(1)-1).*exp(-r*t));hold on;plot(tm2,data2,'g');3：拟合图形如下：（绿色线是预测时的拟合图线，蓝色线是已知数据拟合图线。