微分方程数值解.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 77142005）。实验一 常微分方程初值解法（一）1 实验内容分别用Euler法、改进Euler法、Runge-kutta法求解初值问题u=-2*t*u2u0=12 实验数据与实验结果1）实验结果欧拉法的值与真实的值改进欧拉法的值与真实的值4阶Runge-kutta法的值与真实的值欧拉法的值与真实的值改进欧拉法的值与真实的值4阶Runge-kutta法的值与真实的值从上图可以看出，Euler法计算出的结果与其他方法所得结果相对来说精度太低。4 程序代码清单原式function f=f(t,u)f=-2*t*u*u;欧拉法functiont,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改进欧拉法functiont,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 0x1u0=u1=02 实验数据与实验结果（可用文字描述或贴图的方式进行说明）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) = -h2*sin(i*h); A(i, i+1) = 1; A(i+1, i) = 1;endb(n-1) = -h2*(n-1)*h);Jacobi:functionx,i=Jacobil(A,b,x0,tol,max)n,n=size(A);xold=x0;D=diag(diag(A);E=eye(n);C=E-(DA);d=Db;i=1;while i=max xnew=C*xold+d; if norm(xnew-xold)= to1 x0 = x; x = T*x0 + f; k = k+1; if(k = max) disp(); return; end % disp(k, x); endSor:function x, k = sor(A, b, x0, w, to1, max)if(w =2) erro; return;endD = 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 实验内容分别用向前向后格式求解初值问题ut=2uu2 (0x1,0t1)ux,0=ex (0x1)u0,t=et,u1,t=e1+t (0t1)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.0000x = 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/h2;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); endend向后差分法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);endfor 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层除端点的节点处的值endfor 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);endfunction A,b=matrixback(x0,xf,h,t,N,a)J=(xf-x0)/h;r=a*t/(h2);for n=1:Nb(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; endendendA=diag(1+2*r)*ones(1,J-1); %主对角为1+2r的对角阵for i=1:J-2A(i,i+1)=-r;A(i+1,i)=-r;endfunction 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 endend%用backsub进行回代求解x=backsub(Aug(1:N,1:N),Aug(1:N,N+1);实验四 显格式求解初值解法（四）1 实验内容用显格式求解utt=uxx (0x0)u0,t=u1,t=0 (t0)ux,0=sin4x,ut(x,0)=sin8x (0x1)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*(t2)/h2;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);endfor n=2:N u(1,n)=0; u(J+1,n)=0; for j=2:J u(j,n+1)=r2*(u(j+1,n)+u(j-1,n)+2*(1-r2)*u(j,n)-u(j,n-1); endend