人口问题数据拟合的MATLAB程序.docx
《人口问题数据拟合的MATLAB程序.docx》由会员分享,可在线阅读,更多相关《人口问题数据拟合的MATLAB程序.docx(11页珍藏版)》请在冰豆网上搜索。
![人口问题数据拟合的MATLAB程序.docx](https://file1.bdocx.com/fileroot1/2023-1/30/b9592971-1105-4eb5-9c18-7fd1454b2edb/b9592971-1105-4eb5-9c18-7fd1454b2edb1.gif)
人口问题数据拟合的MATLAB程序
人口问题数据拟合的MATLAB程序
拟合
%拟合数据人口问题
x=[1949195419591964196919741979198419891994];
y=[5.46.06.77.08.19.19.810.311.311.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=A\y'
z=b
(1)+b
(2).*x;
plot(x,y,'r*',x,z),xlabel('x')%用线性回归
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(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*eplot(x,y,'r*',x,z,x,z1,x,z2)%用非线性函数拟合(缺点初值不合适,就得不到解)x=[49545964697479848994];
y=[5.46.06.77.08.19.19.810.311.311.8];
fun=inline('b
(1).*exp(b
(2).*x)','b','x');b0=[20.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=[1949195419591964196919741979198419891994];
y=[5.46.06.77.08.19.19.810.311.311.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=[x19992005]
zz=b
(1)+b
(2).*x1;
y99=[zz(11)-1.96*ezz(11)+1.96*e]y05=[zz(12)-1.96*ezz(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=[2550100200500];
y=[1.93.44.95.66.1];
b0=[23];%参数初值
fun=inline('b
(1).*x./(b
(2)+x)','b','x');%拟合函数
[b,r,j]=nlinfit(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=[1949195419591964196919741979198419891994];
y=[5.466.778.19.19.810.311.111.8];plot(x,y,'.r','MarkerSize',20),xlabel('x')
%解代数方程
[x1,x2]=size(x);
a=[ones(x2,1),x'];
b=a\y'
%线性回归
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).*x)','b','x');
b0=[0.0001,-30],
[deta,r,j]=nlinfit(x',y',fun,b0);
%非线性拟合
xx=[2550100200500];yy=[1.93.44.95.66.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.5^2-w^2)
%用内联函数求零点
fun=inline('51-30*0.5/w-sqrt(0.5^2-w^2)','w');
w0=fzero(fun,0.3)
%用函数文件求零点
functiony=fun(w)
y=51-30*0.5/w-sqrt(0.5^2-w^2);
w0=fzero(@fun,0.3)
%用循环语句搜寻零点
fori=1:
100
w=0.3-0.001*i;
y=30*0.5/w+sqrt(0.5^2-w^2);
ify>=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=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)
t=0:
100;ss=subs(s,'t',t);plot(t,ss),grid
%求数值解
%建立函数文件fun.m:
functiony=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
functionfyy=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)
plot(t,y(:
1),t,y(:
2),'r')
%高阶微分方程边值问题
Z=dsolve('x*D2y-3*Dy=x^2','y
(1)=0','y(5)=0','x')
x=1:
0.01:
5;
s=subs(Z);
plot(x,s)
%lorentz方程组
%周期轨线
functiondot=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,[0200],[1;2;3]);plot3(yy(:
1),yy(:
2),yy(:
3));
%奇怪吸引子
functiondot1=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,[0200],[1;2;3]);plot3(y1(:
1),y1(:
2),y1(:
3),'r');
试卷合理的均衡分配方案的计算机程序(MATLAB)原题:
在大学生数学建模竞赛的评卷工作中,M个评委(M个评委来自不同的学校)要完成N份试卷的打分,竞赛试卷来自K个学校,第j个学校有竞赛试卷lj份,N=?
lj(j,1:
k)为节省人力,每份试卷只要由其中p(p根据回避原则,要求评委不能阅自己学校的试卷,请给出试卷合理的均衡分配方案的数学模型,使各评委阅卷工作量均衡;试卷分配均衡分散。
2。
给出试卷合理的均衡分配方案的计算机程序,要求用MATLAB或C语言编写。
输入参数为p,M,k,n,输出参数为各评委分别阅卷的号码,就下列实例给出问题的答案。
实例:
某省由竞赛试卷568份,32个评委阅卷,40所学校,3程序:
l=[1416102018192510219
191714131211168148615151515151571587
15241515915161517];%我自己随便定的
M=zeros(32,32);
h=1%h评委的代码
j=1
b=0
f=3%f是一份卷子f个人评阅
a=0;
b=l(n);
forn=1:
40
fori=a+1:
b
ifh==n%规避原则
h=h+f
ifh>32
h=h-32;
end
end
M(h,j)=i;
h=h+1;
ifh==33
h=1;
j=j+1
end
end
a=b;
ifn==40
break
end
b=b+l(n+1);
end
M=M%即所求矩阵
求数列2/13/25/38/513/821/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;
fori=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;%赋初值
b=b0;%赋初值
r=r0;%赋初值
n=n0;%赋初值
fori=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;%赋初值
fori=1:
n-1%注意i不能取到10,否则n,10时a(i+1)=a(11).
a(i+1)=r*a(i)+b;%通项公式
end
a%输出a数列的各项值
运行结果
a=
2.330012.980522.033429.728436.269141.828846.554550.571353.985656.8878
张丘建算经百鸡问题及其解(程序)<张丘建算经>百鸡问题及其解(程序)
问题:
鸡翁一,值钱五,鸡母一,值钱三,鸡雏三,值钱一,百钱买百鸡,问鸡翁母雏各几何?
列出方程式:
x+y+z=100
6*x+3*y+(1/3)*z=100
程序1:
只有兩個式子,所以
沒有辦法使用標準的方法去求解。
使用程式設計的方法,是使用暴力法
xyz=[];
forx=0:
100
fory=0:
100
z=100-x-y;
sum=6*x+310*y+(1.0/3.0)*z;
ifsum-100==0
xyz=[xyzxyz];
else
end
end
end
xyz%如果输出的是xyz=[],则表明方程无整数解!
程序2
原方程变形为:
z=100-x-y
y=(100-7*x)/4
编程序如下
x=0;
xyz=[];
whilex<=100
ifmod(100-7*x,4)==0&100-7*x>0
y=(100-7*x)/4;
z=100-x-y;
xyz=[xyzxyz];
end
x=x+1;
end
xyz
-----------------------------------
因为100-7*x>0,即x<14.
所以可将whilex<=100改为whilex<14.这样会提高程序的效率.
对程序2的改进:
程序2'
x=0;
xyz=[];
whilex<14%由z=100-x-y
ifmod(100-7*x,4)==0%由y=(100-7*x)/4
y=(100-7*x)/4;
z=100-x-y;
xyz=[xyzxyz];
end
x=x+1;
end
xyz
计算所有小于1000的Fibonnaci数列%计算所有小于1000的Fibonnaci数
f=[1;1];
i=1;
whilef(i)+f(i+1)<1000f(i+2)=f(i)+f(i+1);i=i+1;
end
f=f',i
%运行结果
f=
Columns1through12
1123581321345589
144
Columns13through16
233377610987i=
15