暑期实习总结报告

一.概要

实习内容包括四部分：微分方程数值解，多元统计分析，最优化以及综合题目。 列表工作概述： 微分方程数值解 学会调用ode45求解常微分方程的初值问题 学会调用bvp4c求解常微分方程的边值问题 学会调用pdepe求解偏微分方程的初边值问题 自己编程部分主要编制了euler法和改进的euler法。 Spss统计分析部分 理解并掌握方差分析的理论和操作 理解相关分析的理论和操作 理解因子分析的理论和操作 理解假设检验的理论和操作 Lingo及最优化 掌握基本操作 会编译一些简单的最优化问题的程序 能看明白一些比较复杂的例子 综合题目

运用数学知识和matlab解决狐兔模型 二.正文

详细总结各部分的内容和所解的问题，包括相关的知识和方法简述，求解的问题，编制的程序，解结果的讨论，对问题进一步的讨论，对自己工作的评价。 第一部分 微分方程数值解部分 1. ode45的调用

1.1调用ode45解决有关传染病的一个例子

写出有关传染病模型的解析表达式如右所示： disi dti(t)s(t)1



i(0)i0运用matlab调用系统函数ode45进行求解，程序如下：

function y =ill(t,x) a=1; b=0.2;

y=[a*x(1)*x(2)-b*x(1);-a*x(1)*x(2)] >> ts=0:50; x0=[0.02,0.98];

[t,x]=ode45('ill',ts,x0);

plot(t,x(:,1),t,x(:,2)),grid 将所得到的各个节点的值作成图形如下：

图形符合模型规律。 1.2 理论解释

ode45方法用来处理非刚性的常微方程初值问题，matlab中调用语句为： [T,Y] = solver(odefun,tspan,y0)

[T,Y] = solver(odefun,tspan,y0,options)

[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options) sol = solver(odefun,[t0 tf],y0...)

其中odefun是这样的一个函数句柄，它表出常微分方程右端表达式

tspan是对时间轴的划分区间 y0表示初始值对应的向量

options一般缺省，表示使用默认值，有特殊情况时使用 sol返回的是一个结构体，通过调用sol.x，sol.y，sol.solver，sol.xe，sol.ye sol.ie来得到sol中计算出的节点值。

T是时间节点对应的向量，Y是解矩阵，TE是起始时间点，YE是初始解，IE是消失的方程的参数i。

1.3.例1：在[1，2]上解微分方程初值问题dy/dx=u/x-（x/u）^2

y(1)=2; 应用ode45函数解得程序是： function dydx=myode45(x,u) dydx=u/x-(x/u)^2; end

[x,u]=ode45(@myode45,[1 2],2); plot(x,u);

得到的数值解对应的图像为：

真解对应的图像为

二者基本一致，可见吻合的很好

bvp4c方法用来处理常微方程边值问题，使用bvp4c时，先要把2阶微分化为1阶，matlab中调用语句为：

sol = bvp4c(odefun,bcfun,solinit)

sol = bvp4c(odefun,bcfun,solinit,options) solinit = bvpinit(x, yinit, params)

其中，odefun是函数句柄，是微分方程dy/dx方程的右端部分。

bcfun是函数句柄，是将边界条件方程全部移到左端，取其左端部分。

solinit是一个对方程初始条件的猜测，可以使用bvpinit函数进行赋值。

options一般缺省，表示使用默认值，有特殊情况时使用sol.x，sol.y，sol.yp与ode45中的意义一致sol.parameters返回带未知参数的常微方程中未知参数的估计值sol.solver存sol得出的解。

solinit函数的赋值语句为 solinit = bvpinit(x, yinit, params)

例如：在[1，2]上解微分方程初值问题y”+|y|=0

y(0)=0;y(4)=-2; 首先将方程降阶得 y2=y1’ y2’=|y1| odefun函数

function dydx = odefun(x,y) dydx = [ y(2)

-abs(y(1))]; bcfun函数

function res = bcfun(ya,yb)

res = [ ya(1)

yb(1) + 2]; 用bvpinit给solinit赋值

solinit = bvpinit(linspace(0,4,5),[1 0]); 调用bvp4c函数

sol = bvp4c(@twoode,@twobc,solinit); 画出对应图形

x = linspace(0,4); y = deval(sol,x); plot(x,y(1,:));

2．偏微方程（pde）数值解

pdepe方法主要解决一维抛物-椭圆方程问题。常用的调用语句是： sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

[sol,tsol,sole,te,ie]=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) 对偏微分问题

c(x,t,u,du/dt)du/dt=x^(-m)*d/dx(x^m*f(x,t,u,du/dx))+s(x,t,u,du/dt) 在t0<=t<=tf,a<=x<=b上处理。 m取值取决于方程的类型，方程是平面型的m=0，方程是圆柱形的m=1，方程是球对称的m=2； 其中pdefun，icfun，bcfun的意义与用法和常微分中的一致。xmesh表示对x区间的划分，

tspan是对t区间的划分。

u2u例如 tt22

u(x,0)sinx u(0,t)0

etu(1,t)0 x程序如下：

表出方程的c，f，s

function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0; 初值函数

function u0 = pdex1ic(x) u0 = sin(pi*x); 边值函数

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0;

pr = pi * exp(-t); qr = 1;

调用pdepe函数 function pdex1

m = 0;

x = linspace(0,1,20); t = linspace(0,2,5);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1);

% A surface plot is often a good way to study a solution. surf(x,t,u)

title('Numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t')

% A solution profile can also be illuminating.

figure

plot(x,u(end,:))

title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)') 得出图像是

以上主要讨论的是方程问题，对方程组问题，只需将各参数令为向量，以同样方式求解即可。 例如:

u12u10.0242F(u1u2) txu22u20.1702F(u1u2) tx其中F(y)exp(5.73y)exp(11.46y)

u2(0,t)0u1(1,t)0 u2(1,t)0x程序如下：

% -------------------------------------------------------------- function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [1; 1];

f = [0.024; 0.17] .* DuDx; y = u(1) - u(2);

F = exp(5.73*y)-exp(-11.47*y); s = [-F; F];

% -------------------------------------------------------------- function u0 = pdex4ic(x); u0 = [1; 0];

% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0];

pr = [ur(1)-1; 0]; qr = [0; 1]; 调用函数是 function pdex4 m = 0;

x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2);

figure

surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t')

figure

surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t') 图像是：

3.微分方程练习题

didtsiidssidt i0i0,s0s0%example11_myEuler N=20; %节点数 h=1/N; %步长

u=zeros(1,N+1); %估计值 x=zeros(1,N+1); %x值向量 u(1)=2; %初始化

for j=1:N+1 %循环计算x向量 x(j)=1+(j-1)*h; end

for j=1:N %Euler法的计算

u(j+1)=u(j)+h*(u(j)/x(j)-(x(j)^2)/(u(j)^2)); end %画图

plot(x,u,'r*') %计算值 hold on

fplot(@(x)x*(8-3*log(x))^(1/3),[1 2]) %真值 hold off

%example12_myModEuler N=500; %节点值 h=1/N; %步长

u=zeros(1,N+1); %计算值 uc=zeros(1,N+1); %真实值 x=zeros(1,N+1); %x向量 %初始化 u(1)=2; uc(1)=2;

u1=zeros(1,6); %迭代向量 for j=1:N+1 %x向量计算 x(j)=1+(j-1)*h; end

for j=1:N %Euler法的估计

u(j+1)=u(j)+h*(u(j)/x(j)-(x(j)^2)/(u(j)^2)); end

for j=1:N %改进的Euler法

u1(1)=u(j)+h*(u(j)/x(j)-(x(j)^2)/(u(j)^2)); for i=1:5 %循环迭代u1处的真值

u1(i+1)=u1(i)+h/2*(u1(i)/x(j+1)-(x(j+1)^2)/(u1(i)^2)+u(j)/x(j)-(x(j)^2)/(u(j)^2)); end

uc(j+1)=u1(6); end %画图 plot(x,uc) hold on

fplot(@(x)x*(8-3*log(x))^(1/3),[1 2]) %真实图形

hold off

h=0.2; tao=0.04;

r=(4*tao)/(h^2)*(pi^2); N=4/h; J=N+1;

u=zeros(N+1,11); for j=1:N+1

u(j,1)=sin(pi*h*(j-1)/4)*(1+2*cos(pi*h*(j-1)/4)); end for i=1:10 for j=2:N if j==2

u(j,i+1)=r*u(j+1,i)+(1-2*r)*u(j,i); elseif j==N+1

u(j,i+1)=(1-2*r)*u(j,i)+r*u(j-1,i); else

u(j,i+1)=r*u(j+1,i)+(1-2*r)*u(j,i)+r*u(j-1,i); end end end

plot(u(2:N,11),'r')

title('parabolic equation!')

% 抛物形方程CN差分 h=0.2; tao=0.04; a=4/pi^2; r=a*tao/h^2; J=4/h+1; N=0.4/tao+1; u=zeros (N,J);

% 边值条件均为0，使用矩阵原始值，不再重新赋值 % 初值条件 for j=1:J

u (1,j)=sin (pi*h*(j-1)/4)*(1+2*cos (pi*h*(j-1)/4)); end

% 构造左边系数矩阵M A=(1+r)*ones(1,J-2); B=(-r/2)*ones(1,J-3);

M=diag(A)+diag(B,1)+diag(B,-1); % 构造右边系数矩阵P A=(1-r)*ones(1,J-2); B=(r/2)*ones(1,J-3);

P=diag(A)+diag(B,1)+diag(B,-1); % 用矩阵向后差分求解

for n=1:N-1

u(n+1,2:J-1)=M^(-1)*P*u (n,2:J-1)'; end

% 三维立体显示

t = linspace(0,0.4,N); x = linspace (0,4,J); figure;

surf(x,t,u(:,:));

title('Numerical solution computed with N*J mesh points.') xlabel('Distance x') ylabel('Time t')

% 画出数值解图与解析值图比较,红色曲线是数值解，黑色为解析解 figure;

fplot ('exp (-1)*sin (pi*x/2)+exp (-0.4/4)*sin (pi*x/4)',[0,4]); hold on;

plot (x,u (N,:),'r')

title ('Solution at t = 0.4') xlabel ('Distance x') ylabel ('u (x,0.4)')

function [y]=myoded(t,x) a=1; b=0.3;

y=[a*x(1)*x(2)-b*x(1);-a*x(2)*x(1)]; end

%---------我是分割线----------------% ts=0:50; x0=[0.02,0.98];

[t,x]=ode45('myoded',ts,x0); plot(t,x(:,1),t,x(:,2)); figure

plot(x(:,2),x(:,1));

2.ode45是解决有关常微分方程初值问题的系统函数，而bvp4c则是用来求解常微分方程边之

第三部分：大题目