人口问题数据拟合的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=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*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');

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 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]=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=[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 11.1 11.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=[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.5^2-w^2)

%用内联函数求零点

fun=inline('51-30*0.5/w-sqrt(0.5^2-w^2)','w');

w0=fzero(fun,0.3)

%用函数文件求零点

function y=fun(w)

y=51-30*0.5/w-sqrt(0.5^2-w^2);

w0=fzero(@fun,0.3)

%用循环语句搜寻零点

for i=1:100

w=0.3-0.001*i;

y=30*0.5/w+sqrt(0.5^2-w^2);

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=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：

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)

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 方程组

%周期轨线

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 dot1=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）为节省人力，每份试卷只要由其中p（p<M<K<<N)个评委进行打分就行，

1。 根据回避原则，要求评委不能阅自己学校的试卷，请给出试卷合理的均衡分配方案的数学模型，使各评委阅卷工作量均衡；试卷分配均衡分散。

2。 给出试卷合理的均衡分配方案的计算机程序，要求用MATLAB或C语言编写。输入参数为p,M,k,n,输出参数为各评委分别阅卷的号码，就下列实例给出问题的答案。

实例：某省由竞赛试卷568份，32个评委阅卷，40所学校，3<lj<30自己设定，给出问题的答案。

程序:

l=[14 16 10 20 18 19 25 10 21 9 19 17 14 13 12 11 16 8 14 8 6 15 15 15 15 15 15 7 15 8 7 15 24 15 15 9 15 16 15 17];%我自己随便定的

M=zeros(32,32);

h=1%h评委的代码

j=1

b=0

f=3 %f是一份卷子f个人评阅

a=0;

b=l(n);

for n=1:40

for i=a+1:b

if h==n %规避原则

h=h+f

if h>32

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

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;%赋初值

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 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*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 x<=100

if mod(100-7*x,4)==0 & 100-7*x>0

y=(100-7*x)/4;

z=100-x-y;

xyz=[xyz x y z];

end

x=x+1;

end

xyz

-----------------------------------

因为 100-7*x>0 ,即x<14.

所以可将 while x<=100 改为while x<14. 这样会提高程序的效率. 对程序2的改进:程序2'

x=0;

xyz=[];

while x<14 %由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

144

Columns 13 through 16

233 377 610 987

i =

15

5 8 13 34 89 21 55

本文来源：https://www.bwwdw.com/article/yfn1.html