微分方程数值解.docx

本科生实验报告

实验课程
微分方程数值解
学院名称
管理科学学院
专业名称
信息与计算科学
学生姓名

学生学号

指导教师
林红霞
实验地点
6C402
实验成绩

二〇 一五 年 十 月 二〇 一五 年 十一 月
填写说明

1、 适用于本科生所有的实验报告（印制实验报告册除外）；
2、 专业填写为专业全称，有专业方向的用小括号标明；
3、 格式要求：
① 用A4纸双面打印（封面双面打印）或在A4大小纸上用蓝黑色水笔书写。
② 打印排版：正文用宋体小四号，1.5倍行距，页边距采取默认形式（上下2.54cm，左右2.54cm，页眉1.5cm，页脚1.75cm）。字符间距为默认值（缩放100%，间距：标准）；页码用小五号字底端居中。
③ 具体要求：
题目（二号黑体居中）；
摘要（“摘要”二字用小二号黑体居中，隔行书写摘要的文字部分，小4号宋体）；
关键词（隔行顶格书写“关键词”三字，提炼3-5个关键词，用分号隔开，小4号黑体)；
正文部分采用三级标题；
第1章 ××(小二号黑体居中，段前0.5行)
1.1 ×××××小三号黑体×××××（段前、段后0.5行）
1.1.1小四号黑体（段前、段后0.5行）
参考文献（黑体小二号居中，段前0.5行），参考文献用五号宋体，参照《参考文献著录规则（GB/T 7714－2005）》。
实验一 常微分方程初值解法（一）

1 实验内容
分别用Euler法、改进Euler法、Runge-kutta法求解初值问题
u·=-2*t*u2u0=1
2 实验数据与实验结果
1）实验结果

欧拉法的值与真实的值

改进欧拉法的值与真实的值

4阶Runge-kutta法的值与真实的值

欧拉法的值与真实的值

改进欧拉法的值与真实的值

4阶Runge-kutta法的值与真实的值

从上图可以看出，Euler法计算出的结果与其他方法所得结果相对来说精度太低。

4 程序代码清单
原式
function f=f(t,u)
f=-2*t*u*u;

欧拉法
function[t,u]=euler(f,t0,u0,tf,h)
n=(tf-t0)/h;
u(1)=u0;
t(1)=t0;
for i=1:n
t(i+1)=t0+i*h;
u(i+1)=u(i)+h*feval(f,t(i),u(i));
end

改进欧拉法
function[t,u]=adeuler(f,t0,u0,tf,h)
% Eluer 方法
% 一阶微分方程的函数: f
% 初始条件： t0, u0
% 取值范围的一个端点：tf
% 区间步长：h(默认值为 0.1)n=fix((tf-t0)/h);
u(1)=u0;
t(1)=t0;
for i=1:n
t(i+1)=t0+i*h;
up=u(i)+h*feval(f,t(i),u(i));
uc=u(i)+h*feval(f,t(i+1),up);
u(i+1)=(up+uc)/2;
end
四阶Runge-kutta法
function [T,Y]=rk4(f,a,b,ya,N)
h=(b-a)/N;
T=a:h:b;
Y(1)=ya;
for j=1:N
k1=h*feval(f,T(j),Y(j));
k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
k4=h*feval(f,T(j)+h,Y(j)+k3);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end

实验二 求解二阶常微分初值问题（二）
1 实验内容
分别用差分法、Jacobi迭代法求解初值问题
-dudx2=sinx 0<x<1u0=u1=0
2 实验数据与实验结果（可用文字描述或贴图的方式进行说明）
1）实验结果

3 程序代码清单（可直接将可运行源代码粘贴在下面的方框中）
Matrix:
function [A,b] = matrix(x0,xf,h)
n = (xf-x0) / h;
A = diag((-2)*ones(1,n-1));
b = zeros(n-1,1);
for i = 1 : n-2
b(i) = -h^2*sin(i*h);
A(i, i+1) = 1;
A(i+1, i) = 1;
end
b(n-1) = -h^2*((n-1)*h);

Jacobi:
function[x,i]=Jacobil(A,b,x0,tol,max)
[n,n]=size(A);
xold=x0;
D=diag(diag(A));
E=eye(n);
C=E-(D\A);
d=D\b;
i=1;
while i<=max
xnew=C*xold+d;
if norm(xnew-xold)<=tol
x=xnew;
disp('jacobiµü´úÊÕÁ²');
return;
else
xold=xnew;
end
%disp([xnew']);
i=i+1;
end
disp('Jacobµü´ú·¨²»ÊÕÁ²');
disp('×î´óµü´ú´ÎÊýºóµÄ½á¹ûΪ£º');
x=xnew;

Gauss:
function [x, k] = Gauss(A, b, x0, to1, max)
D = diag(diag(A));
L = tril(A, -1);
U = triu(A, 1);
T = -(D+L) \ U;
f = (D+L) \ b;
x = T*x0 + f;
k = 1;
while norm(x - x0) >= to1
x0 = x;
x = T*x0 + f;
k = k+1;
if(k >= max)
disp('µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²');
return;
end
% disp([k, x']); ÏÔʾÿ²½½á¹û
end

Sor:
function [x, k] = sor(A, b, x0, w, to1, max)
if(w <= 0 || w >=2)
erro;
return;
end
D = diag(diag(A));
L = tril(A, -1);
U = triu(A, 1);
T = inv(D+L*w)*((1-w)*D-w*U);
f = w*inv(D+L*w)*b;
x = T*x0 + f;
k = 1;
while norm(x - x0) >= to1
x0 = x;
x = T*x0 + f;
k = k+1;
if(k >= max)
disp('µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²');
return;
end
% disp([k, x']); ÏÔʾÿ²½½á¹û
end

实验三 用向前向后差分格式求解（三）

1 实验内容
分别用向前向后格式求解初值问题
∂u∂t=∂2u∂u2 (0≼x≼1,0<t≼1)ux,0=ex (0<x<1)u0,t=et,u1,t=e1+t (0≼t≼1)
2 实验数据与实验结果（可用文字描述或贴图的方式进行说明）
1）实验结果
向前差分结果：
T =
Columns 1 through 10
0 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0.0400 0.0450
Columns 11 through 20
0.0500 0.0550 0.0600 0.0650 0.0700 0.0750 0.0800 0.0850 0.0900 0.0950
Columns 21 through 30
0.1000 0.1050 0.1100 0.1150 0.1200 0.1250 0.1300 0.1350 0.1400 0.1450
Columns 31 through 40
0.1500 0.1550 0.1600 0.1650 0.1700 0.1750 0.1800 0.1850 0.1900 0.1950
Columns 41 through 50
0.2000 0.2050 0.2100 0.2150 0.2200 0.2250 0.2300 0.2350 0.2400 0.2450
Columns 51 through 60
0.2500 0.2550 0.2600 0.2650 0.2700 0.2750 0.2800 0.2850 0.2900 0.2950
Columns 61 through 70
0.3000 0.3050 0.3100 0.3150 0.3200 0.3250 0.3300 0.3350 0.3400 0.3450
Columns 71 through 80
0.3500 0.3550 0.3600 0.3650 0.3700 0.3750 0.3800 0.3850 0.3900 0.3950
Columns 81 through 90
0.4000 0.4050 0.4100 0.4150 0.4200 0.4250 0.4300 0.4350 0.4400 0.4450
Columns 91 through 100
0.4500 0.4550 0.4600 0.4650 0.4700 0.4750 0.4800 0.4850 0.4900 0.4950
Columns 101 through 110
0.5000 0.5050 0.5100 0.5150 0.5200 0.5250 0.5300 0.5350 0.5400 0.5450
Columns 111 through 120
0.5500 0.5550 0.5600 0.5650 0.5700 0.5750 0.5800 0.5850 0.5900 0.5950
Columns 121 through 130
0.6000 0.6050 0.6100 0.6150 0.6200 0.6250 0.6300 0.6350 0.6400 0.6450
Columns 131 through 140
0.6500 0.6550 0.6600 0.6650 0.6700 0.6750 0.6800 0.6850 0.6900 0.6950
Columns 141 through 150
0.7000 0.7050 0.7100 0.7150 0.7200 0.7250 0.7300 0.7350 0.7400 0.7450
Columns 151 through 160
0.7500 0.7550 0.7600 0.7650 0.7700 0.7750 0.7800 0.7850 0.7900 0.7950
Columns 161 through 170
0.8000 0.8050 0.8100 0.8150 0.8200 0.8250 0.8300 0.8350 0.8400 0.8450
Columns 171 through 180
0.8500 0.8550 0.8600 0.8650 0.8700 0.8750 0.8800 0.8850 0.8900 0.8950
Columns 181 through 190
0.9000 0.9050 0.9100 0.9150 0.9200 0.9250 0.9300 0.9350 0.9400 0.9450
Columns 191 through 200
0.9500 0.9550 0.9600 0.9650 0.9700 0.9750 0.9800 0.9850 0.9900 0.9950
Column 201
1.0000
x =
Columns 1 through 10
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000
Column 11
1.0000
向后差分结果：

图3-1 向后差分结果

图3-2 向后差分结果

4 程序代码清单（可直接将可运行源代码粘贴在下面的方框中）
向前差分法
function [T,x,u]=front(a,xa,xb,ta,tb,J,N)
h=(xb-xa)/J;
t=(tb-ta)/N;
r=a*t/h^2;
T(1)=ta;
x(1)=xa;
x(J+1)=xb;
for j=2:J
u(j,1)=exp((j-1)*h); %t0层除端点的节点的值
x(j)=xa+(j-1)*h;
end
u(1,N+1)=exp(N*t); %x0节点处第N层的值
u(J+1,N+1)=exp(1+N*t); %xJ节点处第N层的值
for n=1:N
T(n+1)=ta+n*t;
u(1,n)=exp((n-1)*t); %x0节点处第0层到N-1层的值
u(J+1,n)=exp(1+(n-1)*t); %xJ节点第0层到N-1层的值
for j=2:J
u(j,n+1)=r*u(j+1,n)+(1-2*r)*u(j,n)+r*u(j-1,n);
end
end
向后差分法
function [C,x]=back(A,N,J,h,b)
for j=1:J-1
x(j)=(j-1)*h;
C(j,1)=exp(j*h);
end
for n=2:(N+1)
C(:,n)=uptrbk(A,C(:,n-1)+b(:,n-1));
end

function [C,x]=back(A,N,J,h,b)
%向后差分格式 AU(n+1)=U(n)
for j=1:J-1
x(j)=(j-1)*h;
C(j,1)=exp(j*h); %C的第一列元素,即第0层除端点的节点处的值end
for n=2:(N+1)
C(:,n)=uptrbk(A,C(:,n-1)+b(:,n-1)); %用列主消元法计算
end

function x=backsub(A,b)
%系数矩阵是上三角矩阵求解线性方程组
%从下往求解即解xn,x(n-1),....n=length(b);
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

function [A,b]=matrixback(x0,xf,h,t,N,a)
J=(xf-x0)/h;
r=a*t/(h^2);
for n=1:N
b(1,n)=r*exp(n*t);
b(J-1,n)=r*exp(1+n*t);
for i=2:J-2
for n=1:N
b(i,n)=0;
end
end
end
A=diag((1+2*r)*ones(1,J-1)); %主对角为1+2r的对角阵
for i=1:J-2
A(i,i+1)=-r;
A(i+1,i)=-r;
end

function x=uptrbk(A,b)
[N N]=size(A);
x=zeros(1,N+1);
Aug=[A,b]; % 定义增广矩阵
%将增广矩阵化为上三角，阶梯阵
for p=1:N-1
[Y,j]=max(abs(Aug(p:N,p)));
%增广矩阵中（app,a(p+1)p,....,aNp）绝对值最大记为Y,j是最大值在向量中第j个，则位于矩阵中第j+p-1行
C=Aug(p,:); %C为p行所有元素
Aug(p,:)=Aug(j+p-1,:);%将p行元素与j+p-1行交换；
Aug(j+p-1,:)=C;%记C为增广矩阵j+p-1行元素
if Aug(p,p)==0, % app=0 disp
break;
end
for k=p+1:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);
% Rk行-(akp/app)Rp,将app以下（第p列）化为0
end
end
%用backsub进行回代求解
x=backsub(Aug(1:N,1:N),Aug(1:N,N+1));

实验四 显格式求解初值解法（四）

1 实验内容
用显格式求解
utt=uxx (0<x<1,t>0)u0,t=u1,t=0 (t>0)ux,0=sin4πx,ut(x,0)=sin8πx (0<x<1)
2 实验数据与实验结果（可用文字描述或贴图的方式进行说明）
1） 实验结果

图4-1 (x,u(:,501),'b')图

图4-2 (x,u(:,1001),'b')图

图4-3 (x,u(:,1501),'b')图

4 程序代码清单（可直接将可运行源代码粘贴在下面的方框中）
function [x,u]=bodong(a,xa,xb,ta,tb,h,t)
J=(xb-xa)/h;
N=(tb-ta)/t;
r=a*(t^2)/h^2;
x(1)=xa;
x(J+1)=xa+J*h;
for j=2:J
x(j)=xa+(j-1)*h;
u(j,1)=sin(4*(j-1)*h*pi);
u(j,2)=sin(4*pi*(j-1)*h)+t*sin(8*pi*(j-1)*h);
end
for n=2:N
u(1,n)=0;
u(J+1,n)=0;
for j=2:J
u(j,n+1)=r^2*(u(j+1,n)+u(j-1,n))+2*(1-r^2)*u(j,n)-u(j,n-1);
end
end