微分方程数值解 .pdf
本 科 生 实 验 报 告实 验 课 程微 分 方 程数 值 解学 院 名 称管 理 科 学学 院专 业 名 称信 息 与 计算 科 学学 生学 生 学 号指 导 教 师林 红 霞实 验 地 点6C402 实 验 成 绩二 一五 年 十 月二 一五 年 十一 月精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 20 页填 写 说 明1、 适用于本科生所有的实验报告印制实验报告册除外；2、 专业填写为专业全称，有专业方向的用小括号标明；3、格式要求： 用 A4纸双面打印封面双面打印或在A4大小纸上用蓝黑色水笔书写。 打印排版：正文用宋体小四号，1.5 倍行距，页边距采取默认形式上下 2.54cm，左右 2.54cm，页眉，页脚。字符间距为默认值缩放100% ，间距：标准；页码用小五号字底端居中。 具体要求：题目二号黑体居中；摘要 “摘要”二字用小二号黑体居中，隔行书写摘要的文字部分，小4号宋体；关键词 隔行顶格书写“关键词”三字，提炼3-5 个关键词，用分号隔开，小 4 号黑体) ；正文部分采用三级标题；第 1 章 ( 小二号黑体居中，段前0.5 行) 1.1 小三号黑体小四号黑体参考文献 ，参考文献用五号宋体，参照参考文献著录规则GB/T 77142005 。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 20 页实验一常微分方程初值解法一1 实验内容分别用 Euler 法、改良 Euler 法、Runge-kutta法求解初值问题?= -2 ? ? ?2? (0) = 12 实验数据与实验结果1实验结果欧拉法的值与真实的值改良欧拉法的值与真实的值精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 20 页4 阶Runge-kutta法的值与真实的值欧拉法的值与真实的值改良欧拉法的值与真实的值精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 20 页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;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 20 页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精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 20 页实验二求解二阶常微分初值问题二1 实验内容分别用差分法、 Jacobi迭代法求解初值问题-?2= sin ? 0 ? 1? (0) = ? (1) = 02 实验数据与实验结果可用文字描述或贴图的方式进行说明1实验结果精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 20 页3 程序代码清单可直接将可运行源代码粘贴在下面的方框中精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 20 页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( ?y? ?2?2);return;end% disp(k, x); ?2?1?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( ?y? ?2?2);return;end% disp(k, x); ?2?1?end精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 20 页实验三用向前向后差分格式求解三1 实验内容分别用向前向后格式求解初值问题?=?2?2 (0? ? 1,0 ? ? 1)? (?, 0) = ? (0 ? 1)? (0,? ) = ?,? (1,? ) = ?1+? (0? ? ? 1)2 实验数据与实验结果可用文字描述或贴图的方式进行说明1实验结果向前差分结果：T =Columns 1 through 10Columns 11 through 20Columns 21 through 30Columns 31 through 40Columns 41 through 50Columns 51 through 60Columns 61 through 70Columns 71 through 80Columns 81 through 90Columns 91 through 100Columns 101 through 110Columns 111 through 120Columns 121 through 130Columns 131 through 140Columns 141 through 150Columns 151 through 160Columns 161 through 170精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 20 页Columns 171 through 180Columns 181 through 190Columns 191 through 200Column 201x =Columns 1 through 10Column 11向后差分结果：图 3-1 向后差分结果图 3-2 向后差分结果精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 20 页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); 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精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 20 页 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);% 用列主消元法计算endfunction 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-2for 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行交换；精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 20 页 Aug(j+p-1,:)=C;%记C为增广矩阵 j+p-1行元素 if Aug(p,p)=0, % app=0 dispbreak;endfor 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列化为 0endend%用backsub进行回代求解x=backsub(Aug(1:N,1:N),Aug(1:N,N+1);实验四显格式求解初值解法四精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 20 页1 实验内容用显格式求解?= ? (0 ? 0)? (0,? ) = ? (1,? ) = 0 (? 0)? (?, 0) = ?4?,?(?,0) = ?8? (0 ? 1)2 实验数据与实验结果可用文字描述或贴图的方式进行说明1）实验结果精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 20 页精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 20 页图 4-1 (x,u( :,501),b)图图 4-2 (x,u(:,1001),b)图精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 20 页图 4-3 (x,u(:,1501),b)图4 程序代码清单可直接将可运行源代码粘贴在下面的方框中精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 20 页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精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 20 页
