人口问题数据拟合的MATLAB程序.docx

1、人口问题数据拟合的MATLAB程序人口问题数据拟合的MATLAB程序拟合 %拟合数据 人口问题 x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994; y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8; % 1 线性模型 %用一阶多项式 b=polyfit(x,y,1) z=b(2)+b(1).*x; plot(x,y,r*,x,z),xlabel(x) %用矩阵运算 A=ones(size(x), x; b=Ay z=b(1)+b(2).*x; plot(x,y,r*,x,z),xlabel(x) %用

2、线性回归 A=ones(size(x), x; b,c,r,j,R =regress(y,A) % b 回归系数 c 回归系数的置信区间 r 残差 j 拟合数据的置信区间 R 相关系数 F值、p值 z=b(1)+b(2).*x; z1=z+j(:,1); z2=z+j(:,2); plot(x,y,r*,x,z,x,z1,x,z2),xlabel(x) e=sqrt(sum(y-z).2)/8) zz1=z-1.96*e; zz2=z+1.96*e; plot(x,y,r*,x,z,x,zz1,x,zz2) % 2 非线性模型 y=b(2)exp(b(1)x) %转化为线性函数 A=ones(

3、size(x), x; y1=log(y); b1,r,j,R=regress(y1,A) b=exp(b1(1) b1(2) z=b(1).*exp(b(2).*x); e=sqrt(sum(y-z).2)/8) z1=z-1.96*e; z2=z+1.96*e plot(x,y,r*,x,z,x,z1,x,z2) %用非线性函数拟合(缺点 初值不合适，就得不到解) x=49 54 59 64 69 74 79 84 89 94; y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8; fun=inline(b(1).*exp(b(2).*x),b,x);

4、 b0=2 0.01; b,r,j=nlinfit(x,y,fun,b0) z=b(1).*exp(b(2).*x); plot(x,y,r*,x,z) nlintool(x,y,fun,b0) %拟合曲线图。 %预报 x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994; y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8; A=ones(size(x), x; b,c,r,j,R =regress(y,A) z=b(1)+b(2).*x; e=sqrt(sum(z-y).2)/8) x1=x 1999

5、2005 zz=b(1)+b(2).*x1; y99=zz(11)-1.96*e zz(11)+1.96*e y05=zz(12)-1.96*e zz(12)+1.96*e z1=z+j(:,1) y99(1) y05(1); z2=z+j(:,2) y99(2) y05(2); plot(x,y,r*,x1,zz,x1,z1,x1,z2),xlabel(x) %人员疏散问题 x=25 50 100 200 500; y=1.9 3.4 4.9 5.6 6.1; b0=2 3; %参数初值 fun=inline( b(1).*x./(b(2)+x),b,x); %拟合函数 b,r,j=nlin

6、fit(x,y,fun,b0) %拟合函数的系数、残差 z=b(1).*x./(b(2)+x); z1=z+j(:,1); z2=z+j(:,2); plot(x,y,*r,x,z,x,z1,x,z2) e=sqrt(sum(y-z).2)/3) zz1=z-1.96*e; zz2=z+1.96*e; plot(x,y,*r,x,z,r,x,zz1,x,zz2) nlintool(x,y,fun,b0) %拟合曲线图。 clf x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994; y=5.4 6 6.7 7 8.1 9.1 9.8 10.3

7、11.1 11.8; plot(x,y,.r,MarkerSize,20),xlabel(x) %解代数方程 x1,x2=size(x); a=ones(x2,1),x; b=ay %线性回归 b=regress(y,a) y1=b(1)+b(2).*x n=1;%一阶多项式模型 P=polyfit(x,y,n) y2=polyval(P,x) %指数模型转化为线性模型 yl=log(y); Pl=polyfit(x,yl,1) y3=exp(Pl(2).*exp(Pl(1).*x) b2=exp(Pl(2) %指数模型非线性拟合(失败) fun=inline(b(1).*exp(b(2).*

8、x),b,x); b0=0.0001,-30, deta,r,j=nlinfit(x,y,fun,b0); %非线性拟合 xx=25 50 100 200 500; yy=1.9 3.4 4.9 5.6 6.1; beta0=2,3; fun1=inline(b(1).*xx./(b(2)+xx),b,xx); beta,r,j=nlinfit(xx,yy,fun1,beta0) plot(x,y1,x,y2,x,y3,x,y,.r,MarkerSize,20),xlabel(x) 解方程零点的MATLAB语言 零点 % 解方程 51=30*0.5/w+sqrt(0.52-w2) %用内联函数

9、求零点 fun=inline(51-30*0.5/w-sqrt(0.52-w2),w); w0=fzero(fun,0.3) %用函数文件求零点 function y=fun(w) y=51-30*0.5/w-sqrt(0.52-w2); w0=fzero(fun,0.3) %用循环语句搜寻零点 for i=1:100 w=0.3-0.001*i; y=30*0.5/w+sqrt(0.52-w2); if y =51 w0=w, break; end; end; w0 %一元函数直接赋值做图 w=0.2:0.001:0.4; y=51-30.*0.5./w-sqrt(0.5.2-w.2); x

10、=zeros(size(w); plot(w,y,w,x),grid %一元函数用内联函数赋值做图 y=inline(-51+30*0.5/w+sqrt(0.5*2-w*2),w); w=0.1:0.01:0.5; y_char=vectorize(y); Y=feval(y_char,w); plot(w,Y) 求解微分方程的MATLAB语言程序 微分方程 % 一阶微分方程 % 求表达式(符号运算) S=dsolve(Dx=r*x*(1-x/k) S1=dsolve(Dx=r*x*(1-x/k),x(0)=0.01) % 转换为函数值 r=0.3;k=8;C1=0.01; s=subs(S)

11、 t=0:100; ss=subs(s,t,t); plot(t,ss),grid % 求数值解 % 建立函数文件fun.m: function y=fun(x,t) y=0.3*x*(1-x/8); %调用ode23(二、三阶龙格库塔法)解微分方程 t,x=ode23(fun,0, 20,0.01) plot(t,x) %高阶微分方程初值问题 %求数值解要化为一阶方程组 %建立M文件 funs. M function fyy=funs(t,y) fyy=y(1).*(1-y(2).2)-y(2);y(1); t,y=ode23(funs,0,25, 0.25,0 ; plot(t,y) pl

12、ot(t,y(:,1),t,y(:,2),r) %高阶微分方程边值问题 Z=dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x) x=1:0.01:5; s=subs(Z); plot(x,s) % lorentz 方程组 %周期轨线 function dot=odone(t,y) dot=10*(y(2)-y(1); 28*y(1)-y(2)-y(1)*y(3); y(1)*y(2)-8/3*y(3); tt,yy=ode45(odone,0 200,1; 2; 3); plot3(yy(:,1),yy(:,2),yy(:,3); %奇怪吸引子 function dot

13、1=odoneq(t,y) dot1=10*(y(2)-y(1)-3*(28*y(1)-y(2)-y(1)*y(3); 28*y(1)-y(2)-y(1)*y(3); y(1)*y(2)-8/3*y(3); t1,y1=ode45(odoneq,0 200,1; 2; 3); plot3(y1(:,1),y1(:,2),y1(:,3),r); 试卷合理的均衡分配方案的计算机程序(MATLAB) 原题: 在大学生数学建模竞赛的评卷工作中，M个评委(M个评委来自不同的学校)要完成N份试卷的打分，竞赛试卷来自K个学校，第j个学校有竞赛试卷lj份，N=?lj(j,1:k)为节省人力，每份试卷只要由其中

14、p(pMKN)个评委进行打分就行， 1。 根据回避原则，要求评委不能阅自己学校的试卷，请给出试卷合理的均衡分配方案的数学模型，使各评委阅卷工作量均衡;试卷分配均衡分散。 2。 给出试卷合理的均衡分配方案的计算机程序，要求用MATLAB或C语言编写。输入参数为p,M,k,n,输出参数为各评委分别阅卷的号码，就下列实例给出问题的答案。 实例:某省由竞赛试卷568份，32个评委阅卷，40所学校，3lj32 h=h-32; end end M(h,j)=i; h=h+1; if h=33 h=1; j=j+1 end end a=b; if n=40 break end b=b+l(n+1); end

15、 M=M %即所求矩阵 求数列2/1 3/2 5/3 8/5 13 /8 21/13 。的前n(,100)项和 求数列2/1， 3/2 ，5/3 ，8/5 ，13 /8 ，21/13 。的前n(,100)项和 a(1)=2;a(2)=3; b(1)=1;b(2)=2; n=100; sum=0; for i=1:n a(i+2)=a(i)+a(i+1); b(i+2)=b(i)+b(i+1); sum=sum+a(i)/b(i); end sum 结果 sum = 162.1030 差分方程的matlab解法 差分方程的一般形式为:a(n+1),r*a(n)，b 计算程序: a(1)=a0;%

16、赋初值 b=b0;%赋初值 r=r0;%赋初值 n=n0;%赋初值 for i=1:n-1 a(i+1)=r*a(i)+b; %通项公式 end a %输出a数列的各项值 实例: 比如要计算差分方程 a(n+1),0.85*a(n)，11，a(1),2.33的前10项， 可写入下列代码: a(1)=2.33;%赋初值 b=11;%赋初值 r=0.85;%赋初值 n=10;%赋初值 for i=1:n-1 %注意i不能取到10，否则n,10时a(i+1)=a(11). a(i+1)=r*a(i)+b; %通项公式 end a %输出a数列的各项值 运行结果 a = 2.3300 12.9805

17、22.0334 29.7284 36.2691 41.8288 46.5545 50.5713 53.9856 56.8878 张丘建算经 百鸡问题及其解(程序) 百鸡问题及其解(程序) 问题: 鸡翁一,值钱五,鸡母一,值钱三,鸡雏三,值钱一,百钱买百鸡,问鸡翁母雏各几何? 列出方程式 : x + y + z = 100 6*x + 3*y + (1/3)*z= 100 程序1: 只有兩個式子，所以 沒有辦法使用標準的方法去求解。 使用 程式設計的方法，是 使用暴力法 xyz=; for x=0 :100 for y=0 :100 z= 100 - x - y; sum= 6*x + 310*

18、y + (1.0/3.0)*z; if sum-100=0 xyz=xyz x y z; else end end end xyz %如果输出的是xyz=,则表明方程无整数解! 程序2 原方程变形为: z = 100-x-y y=(100-7*x)/4 编程序如下 x=0; xyz=; while x0 y=(100-7*x)/4; z=100-x-y; xyz=xyz x y z; end x=x+1; end xyz - 因为 100-7*x0 ,即x14. 所以可将 while x=100 改为while x14. 这样会提高程序的效率. 对程序2的改进:程序2 x=0; xyz=; w

19、hile x14 %由z = 100-x-y if mod(100-7*x,4)=0 %由y=(100-7*x)/4 y=(100-7*x)/4; z=100-x-y; xyz=xyz x y z; end x=x+1; end xyz 计算所有小于1000的Fibonnaci数列 % 计算所有小于1000的Fibonnaci数 f=1;1; i=1; while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1; end f=f,i %运行结果 f = Columns 1 through 12 1 1 2 3 5 8 13 21 34 55 89 144 Columns 13 through 16 233 377 610 987 i = 15