当前位置:首页 >> 数学 >>

偏微分方程数值解法题解


偏微分方程数值解法(带程序)
? ?u ? 2 u ? 2 , x ? (0,1), t ? 0 ? ? t ?x ? ? 1 ? 2 x, x ? (0, ] ? ? 例 1 求解初边值问题 ? ? 2 要求采用树脂格式 u ( x,0) ? ? ? ?2(1 ? x), x ? [ 1 ,1) ? ? 2 ? ?

u (0, t ) ? u (1, t ) ? 0, t ? 0
?1 n n n un ? un j j ? ? (u j ?1 ? 2u j ? u j ?1 ) , ? ?

?t ,完成下列计算: (?x)2

(1) 取 ?x ? 0.1, ? ? 0.1, 分别计算 t ? 0.01,0.02,0.1, 时刻的数值解。 (2) 取 ?x ? 0.1, ? ? 0.5, 分别计算 t ? 0.01,0.02,0.1, 时刻的数值解。 (3) 取 ?x ? 0.1, ? ? 1.0, 分别计算 t ? 0.01,0.02,0.1, 时刻的数值解。 并与解析解 u( x, t ) ?

8

?

2

sin( )sin(n? x)e( ? n ? t ) 进行比较。 ? 2 2 n n ?1
2 2

?

1

n?

解:程序 function A=zhongxinchafen(x,y,la) U=zeros(length(x),length(y)); for i=1:size(x,2) if x(i)>0&x(i)<=0.5 U(i,1)=2*x(i); elseif x(i)>0.5&x(i)<1 U(i,1)=2*(1-x(i)); end end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)+la*(U(i+2,j)-2*U(i+1,j)+U(i,j)); end end A=U(:,size(U,2)) function u=jiexijie1(x,t) for i=1:size(x,2) k=3;

a1=(1/(1^2)*sin(1*pi/2)*sin(1*pi*x(i))*exp(-1^2*pi^2*t)); a2=a1+(1/(2^2)*sin(2*pi/2)*sin(2*pi*x(i))*exp(-2^2*pi^2*t)); while abs(a2-a1)>0.00001 a1=a2; a2=a1+(1/(k^2)*sin(k*pi/2)*sin(k*pi*x(i))*exp(-k^2*pi^2*t)); k=k+1; end u(i)=8/(pi^2)*a2; end clc; %第 1 题第 1 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.001:t1];y2=[0:0.001:t2];y3=[0:0.001:t3];la=0.1; subplot(131) A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5);line(x,u3,'color','b','linewidth',1); title('例 1(1)'); subplot(132); line(x,u1,'color','b','linewidth',1);line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); clc; %第 1 题第 2 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.005:t1];y2=[0:0.005:t2];y3=[0:0.005:t3];la=0.5; subplot(131); A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3)
-1-

line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例 1(2)'); subplot(132); line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); clc; %第 1 题第 3 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.01:t1];y2=[0:0.01:t2];y3=[0:0.01:t3];la=1.0; subplot(131); A1=zhongxinchafen(x,y1,la);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=zhongxinchafen(x,y2,la);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); A3=zhongxinchafen(x,y3,la);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例 1(3)'); subplot(132); line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); 运行结果:
表 1:取 ?x ? 0.1, ? ? 0.1, t ? 0.01,0.02,0.1, 时刻的解析解与数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984

0 0.1996 0.3968 0.5822 0.7281 0.7867 0.7281

t ? 0.02

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328
-2-

0 0.1939 0.3781 0.5373 0.6487 0.6891 0.6487

t ? 0.1

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873

0 0.0944 0.1796 0.2472 0.2907 0.3056 0.2907

0.5941 0.4317 0.2269 0.0000

0.5822 0.3968 0.1996 0

0.5383 0.3911 0.2056 0.0000

0.5373 0.3781 0.1939 0

0.2444 0.1776 0.0934 0.0000

0.2472 0.1796 0.0944 0

表 2:取 ?x ? 0.1, ? ? 0.5, t ? 0.01,0.02,0.1, 时刻的解析解与数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984 0.5941 0.4317 0.2269 0.0000

0
0.2000 0.4000 0.6000 0.7000 0.8000 0.7000 0.6000 0.4000 0.2000

t ? 0.02

0

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328 0.5383 0.3911 0.2056 0.0000

0
0.2000 0.3750 0.5500 0.6250 0.7000 0.6250 0.5500 0.3750 0.2000

t ? 0.1

0

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873 0.2444 0.1776 0.0934 0.0000

0
0.0949 0.1717 0.2484 0.2778 0.3071 0.2778 0.2484 0.1717 0.0949

0

表 3:取 ?x ? 0.1, ? ? 1.0, t ? 0.01,0.02,0.1, 时刻的解析解与数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984 0.5941 0.4317 0.2269 0.0000

0 0.2000 0.4000 0.6000 0.8000 0.6000 0.8000 0.6000 0.4000 0.2000 0

t ? 0.02

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328 0.5383 0.3911 0.2056 0.0000

0 0.2000 0.4000 0.6000 0.4000 1.0000 0.4000 0.6000 0.4000 0.2000 0

t ? 0.1

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873 0.2444 0.1776 0.0934 0.0000

0 220.200 -453.600 684.600 -861.200 929.000 -861.200 684.600 -453.600 220.200 0

-3-

图 1:取 ?x ? 0.1, ? ? 0.1, t ? 0.01,0.02,0.1, 时刻的解析解与数值解

图2:取 ?x ? 0.1, ? ? 0.5, t ? 0.01,0.02,0.1, 时刻的解析解与数值解

-4-

图3:取 ?x ? 0.1, ? ? 1.0, t ? 0.01,0.02,0.1, 时刻的解析解与数值解

例 2 用 Crank-Nicolson 格式完成例 1 的所有任务。 解: Crank ? Nicolson 格式
?1 un ? un j j

?
? a

?

?1 ?1 ?1 ? (u n ? 2u n ? un ) ? (u n ? 2u n ? un )? ? 0 j ?1 j j ?1 j ?1 j j ?1 ? ? 2h 2

a

将第 n 层和第 n ? 1 层放到等式两端,得

1 a n ?1 a n ?1 a n 1 a n a n n ?1 u ? ( ? ) u ? u ? u ? ( ? ) u ? u j ?1 j ?1 j j ?1 j ?1 j 2 2 2 2 ? h2 ? h2 2h 2h 2h 2h
? , ? ? ? * h2 ( ? 在程序中用 la 表示) , ? 为时间步长, h 为空间步长(在 2 h

??

程序中用 dt , dx 来表示) 化简:左右两边同时乘以 2a? 得:
?1 n ?1 n ?1 n n n ?a?u n ? (2 ? 2a? )u j ? a?u j ?1 ? a?u j ?1 ? (2 ? 2a? )u j ? a?u j ?1 j ?1

程序: function A=granknicolson(x,y,la,a) U=zeros(size(x,2),size(y,2)); for i=1:size(x,2) if x(i)>0&x(i)<=0.5 U(i,1)=2*x(i); elseif x(i)>0.5&x(i)<1 U(i,1)=2*(1-x(i)); end
-5-

end x=[0:0.1:1]; V1=zeros(size(x,2)-2); V2=zeros(size(x,2)-2); for i=1:size(x,2)-2 for j=1:size(x,2)-2 if i==j V1(i,j)=2+2*a*la; V2(i,j)=2-2*a*la; elseif i==j+1|j==i+1 V1(i,j)=-a*la; V2(i,j)=a*la; end end end for j=2:size(y,2) U(2:size(x,2)-1,j)=inv(V1)*V2*U(2:size(x,2)-1,j-1); end A=U(:,size(U,2)); function u=jiexijie1(x,t) for i=1:size(x,2) k=3; a1=(1/(1^2)*sin(1*pi/2)*sin(1*pi*x(i))*exp(-1^2*pi^2*t)); a2=a1+(1/(2^2)*sin(2*pi/2)*sin(2*pi*x(i))*exp(-2^2*pi^2*t)); while abs(a2-a1)>0.00001 a1=a2; a2=a1+(1/(k^2)*sin(k*pi/2)*sin(k*pi*x(i))*exp(-k^2*pi^2*t)); k=k+1; end u(i)=8/(pi^2)*a2; end clc; %第 2 题第 1 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.001:t1];y2=[0:0.001:t2];y3=[0:0.001:t3];la=0.1;a=1; subplot(131); A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1);
-6-

A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例 2(1)'); subplot(132); line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); clc; %第 2 题第 2 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.005:t1];y2=[0:0.005:t2];y3=[0:0.005:t3];la=0.5;a=1; subplot(131); A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2) line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例 2(2)'); subplot(132); line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); clc; %第 2 题第 3 问 clear; t1=0.01;t2=0.02;t3=0.1;x=[0:0.1:1]; y1=[0:0.01:t1];y2=[0:0.01:t2];y3=[0:0.01:t3];la=1.0;a=1; subplot(131); A1=granknicolson(x,y1,la,a);u1=jiexijie1(x,t1) line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); A2=granknicolson(x,y2,la,a);u2=jiexijie1(x,t2)
-7-

line(x,A2,'color','r','linestyle',':','linewidth',1.5);line(x,u2,'color','b','linewidth',1); A3=granknicolson(x,y3,la,a);u3=jiexijie1(x,t3) line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1);title('例 2(3)'); subplot(132); line(x,u1,'color','b','linewidth',1); line(x,u2,'color','b','linewidth',1); line(x,u3,'color','b','linewidth',1);title('解析解'); subplot(133); line(x,A1,'color','r','linestyle',':','linewidth',1.5); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,A3,'color','r','linestyle',':','linewidth',1.5);title('数值解'); 运行结果:
表 4:取 ?x ? 0.1, ? ? 0.1, t ? 0.01,0.02,0.1, 时刻的解析解与 Crank-Nicolson 格式数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984 0.5941 0.4317 0.2269 0.0000

0 0.1993 0.3959 0.5810 0.7288 0.7904 0.7288 0.5810 0.3959 0.1993 0

t ? 0.02

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328 0.5383 0.3911 0.2056 0.0000

0 0.1933 0.3773 0.5371 0.6500 0.6914 0.6500 0.5371 0.3773 0.1933 0

t ? 0.1

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873 0.2444 0.1776 0.0934 0.0000

0 0.0949 0.1805 0.2484 0.2921 0.3071 0.2921 0.2484 0.1805 0.0949 0

表 5:取 ?x ? 0.1, ? ? 0.5, t ? 0.01,0.02,0.1, 时刻的解析解与 Crank-Nicolson 格式数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984 0.5941 0.4317 0.2269

0
0.1992 0.3959 0.5820 0.7293 0.7879 0.7293 0.5820 0.3959 0.1992

t ? 0.02

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328 0.5383 0.3911 0.2056
-8-

0
0.1934 0.3776 0.5375 0.6497 0.6906 0.6497 0.5375 0.3776 0.1934

t ? 0.1

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873 0.2444 0.1776 0.0934

0
0.0949 0.1804 0.2484 0.2920 0.3070 0.2920 0.2484 0.1804 0.0949

0.0000

0

0.0000

0

0.0000

0

表 6:取 ?x ? 0.1, ? ? 1.0, t ? 0.01,0.02,0.1, 时刻的解析解与 Crank-Nicolson 格式数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解 时间(t) 解析解 数值解

0 0.2269 0.4317 0.5941 t ? 0.01 0.6984 0.7344 0.6984 0.5941 0.4317 0.2269 0.0000

0 0.1989 0.3956 0.5834 0.7381 0.7691 0.7381 0.5834 0.3956 0.1989 0

t ? 0.02

0 0.2056 0.3911 0.5383 0.6328 0.6654 0.6328 0.5383 0.3911 0.2056 0.0000

0 0.1936 0.3789 0.5397 0.6461 0.6921 0.6461 0.5397 0.3789 0.1936 0

t ? 0.1

0 0.0934 0.1776 0.2444 0.2873 0.3021 0.2873 0.2444 0.1776 0.0934 0.0000

0 0.0948 0.1803 0.2482 0.2918 0.3069 0.2918 0.2482 0.1803 0.0948 0

图4:取 ?x ? 0.1, ? ? 0.1, t ? 0.01,0.02,0.1, 时刻的解析解与Crank-Nicolson格式数值解

-9-

图5:取 ?x ? 0.1, ? ? 0.5, t ? 0.01,0.02,0.1, 时刻的解析解与Crank-Nicolson格式数值解

图6:取 ?x ? 0.1, ? ? 1.0, t ? 0.01,0.02,0.1, 时刻的解析解与Crank-Nicolson格式数值解

?ut ? u xx , x ? (0,1), t ? 0 ? u ( x, 0) ? 1, x ? (0,1) ? ? ?u ? u x ?0 , t ? 0 例 3 数值求解初边值问题 ? ? x x ? 0 ? ? ?u ? ? u x ?1 , t ? 0 ? ? ?x x ?1

要求边界条件离散采用 (1) 向前(后)一阶差商代替微商。 (2) 中心差商(二阶)代替微商。
- 10 -

?1 n n n 数值格式仍采用 un ? un j j ? ? (u j ?1 ? 2u j ? u j ?1 ) , ? ?

?t , (?x)2

这里取 ?x ? 0.1, ? ? 0.25, 计算 t ? 0.005,0.01,0.015,0.02,0.10,0.20,0.25,0.50,1.00
? ? sec ? n ?4?n2t 1 ? 注:本例的解析解为 u ( x, t ) ? 4? ? e cos(2 ? ( x ? )) ?, x ? (0,1) n 2 2 ? n ?1 ? 3 ? 4? n

其中 ? n 是 ? tan ? ?

1 的正根。 2

解:程序: function A=yijiechashang(x,y,la) U=zeros(length(x),length(y)); for i=2:length(x)-1 U(i,1)=1; end U(1,1)=U(2,1)/1.1; U(length(x),1)=U(length(x)-1,1)/1.1; for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)+la*(U(i+2,j)-2*U(i+1,j)+U(i,j)); U(1,j+1)=U(2,j+1)/1.1; U(length(x),j+1)=U(length(x)-1,j+1)/1.1; end end A=U(:,size(U,2)); function u=jiexijie2(x,t,d) r=zhenggen(d); for i=1:length(x) k=3; a1=exp(-4*r(1)^2*t)*cos(2*r(1)*(x(i)-0.5))/((3+4*r(1)^2)*cos(r(1))); a2=a1+exp(-4*r(2)^2*t)*cos(2*r(2)*(x(i)-0.5))/((3+4*r(2)^2)*cos(r(2))); while abs(a2-a1)>0.00001 a1=a2; a2=a1+exp(-4*r(k)^2*t)*cos(2*r(k)*(x(i)-0.5))/((3+4*r(k)^2)*cos(r(k))); k=k+1; end u(i)=4*a2; end clc; %第 3 题第 1 问 clear; t1=0.005;t2=0.01;t3=0.015;t4=0.02;t5=0.10;t6=0.20;t7=0.25;t8=0.5;t9=1.0;
- 11 -

la=0.25;x=[0:0.1:1]; y1=[0:0.0025:t1];y2=[0:0.0025:t2];y3=[0:0.0025:t3];y4=[0:0.0025:t4]; y5=[0:0.0025:t5];y6=[0:0.0025:t6];y7=[0:0.0025:t7];y8=[0:0.0025:t8]; y9=[0:0.0025:t9];d=0:0.000005:100; A1=yijiechashang(x,y1,la);u1=jiexijie2(x,t1,d); A2=yijiechashang(x,y2,la);u2=jiexijie2(x,t2,d); A3=yijiechashang(x,y3,la);u3=jiexijie2(x,t3,d); A4=yijiechashang(x,y4,la);u4=jiexijie2(x,t4,d); A5=yijiechashang(x,y5,la);u5=jiexijie2(x,t5,d); A6=yijiechashang(x,y6,la);u6=jiexijie2(x,t6,d); A7=yijiechashang(x,y7,la);u7=jiexijie2(x,t7,d); A8=yijiechashang(x,y8,la);u8=jiexijie2(x,t8,d); A9=yijiechashang(x,y9,la);u9=jiexijie2(x,t9,d); subplot(121); line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1); line(x,A4,'color','r','linestyle',':','linewidth',1.5); line(x,u4,'color','b','linewidth',1);title('例 3(1)');legend('数值解','解析解'); subplot(122); line(x,A5,'color','r','linestyle',':','linewidth',1.5); line(x,u5,'color','b','linewidth',1); line(x,A6,'color','r','linestyle',':','linewidth',1.5); line(x,u6,'color','b','linewidth',1); line(x,A7,'color','r','linestyle',':','linewidth',1.5); line(x,u7,'color','b','linewidth',1); line(x,A8,'color','r','linestyle',':','linewidth',1.5); line(x,u8,'color','b','linewidth',1); line(x,A9,'color','r','linestyle',':','linewidth',1.5); line(x,u9,'color','b','linewidth',1);legend('数值解','解析解'); clc; %第 3 题第 2 问 clear; t1=0.005;t2=0.01;t3=0.015;t4=0.02;t5=0.10;t6=0.20;t7=0.25;t8=0.5;t9=1.0; la=0.25;x=[0:0.1:1]; y1=[0:0.0025:t1];y2=[0:0.0025:t2];y3=[0:0.0025:t3];y4=[0:0.0025:t4]; y5=[0:0.0025:t5];y6=[0:0.0025:t6];y7=[0:0.0025:t7];y8=[0:0.0025:t8]; y9=[0:0.0025:t9];d=0:0.000005:100; A1=zhongxinchafen2(x,y1,la);u1=jiexijie2(x,t1,d); A2=zhongxinchafen2(x,y2,la);u2=jiexijie2(x,t2,d); A3=zhongxinchafen2(x,y3,la);u3=jiexijie2(x,t3,d);
- 12 -

A4=zhongxinchafen2(x,y4,la);u4=jiexijie2(x,t4,d); A5=zhongxinchafen2(x,y5,la);u5=jiexijie2(x,t5,d); A6=zhongxinchafen2(x,y6,la);u6=jiexijie2(x,t6,d); A7=zhongxinchafen2(x,y7,la);u7=jiexijie2(x,t7,d); A8=zhongxinchafen2(x,y8,la);u8=jiexijie2(x,t8,d); A9=zhongxinchafen2(x,y9,la);u9=jiexijie2(x,t9,d); subplot(121); line(x,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x,u1,'color','b','linewidth',1); line(x,A2,'color','r','linestyle',':','linewidth',1.5); line(x,u2,'color','b','linewidth',1); line(x,A3,'color','r','linestyle',':','linewidth',1.5); line(x,u3,'color','b','linewidth',1); line(x,A4,'color','r','linestyle',':','linewidth',1.5); line(x,u4,'color','b','linewidth',1);title('例 3(2)');legend('数值解','解析解'); subplot(122); line(x,A5,'color','r','linestyle',':','linewidth',1.5); line(x,u5,'color','b','linewidth',1); line(x,A6,'color','r','linestyle',':','linewidth',1.5); line(x,u6,'color','b','linewidth',1); line(x,A7,'color','r','linestyle',':','linewidth',1.5); line(x,u7,'color','b','linewidth',1); line(x,A8,'color','r','linestyle',':','linewidth',1.5); line(x,u8,'color','b','linewidth',1); line(x,A9,'color','r','linestyle',':','linewidth',1.5); line(x,u9,'color','b','linewidth',1);legend('数值解','解析解'); 运行结果:
表 7:取 ?x ? 0.1, ? ? 0.25, 用向前(后)一阶差商代替微商计算

t ? 0.005,0.01,0.015,0.02,0.10,0.20,0.25,0.50,1.00 时刻的解析解与数值解
时间(t) 解析解 数值解 时间 解析解 数值解 时间(t) 解析解 数值解

0.005

0.9250 0.9841 0.9984 0.9999 1.0000 1.0000 1.0000 0.9999 0.9984 0.9841 0.9250

0.8734 0.9607 0.9943 1.0000 1.0000 1.0000 1.0000 1.0000 0.9943 0.9607 0.8734

0.01

0.8965 0.9627 0.9905 0.9984 0.9998 1.0000 0.9998 0.9984 0.9905 0.9627 0.8965
- 13 -

0.8507 0.9358 0.9801 0.9961 0.9996 1.0000 0.9996 0.9961 0.9801 0.9358 0.8507

0.015

0.8755 0.9444 0.9802 0.9945 0.9988 0.9996 0.9988 0.9945 0.9802 0.9444 0.8755

0.8331 0.9164 0.9662 0.9895 0.9976 0.9993 0.9976 0.9895 0.9662 0.9164 0.8331

0.02

0.25

0.8585 0.9286 0.9695 0.9891 0.9967 0.9985 0.9967 0.9891 0.9695 0.9286 0.8585 0.5546 0.6052 0.6454 0.6747 0.6924 0.6984 0.6924 0.6747 0.6454 0.6052 0.5546

0.8184 0.9003 0.9532 0.9817 0.9941 0.9973 0.9941 0.9817 0.9532 0.9003 0.8184 0.5206 0.5727 0.6142 0.6444 0.6628 0.6689 0.6628 0.6444 0.6142 0.5727 0.5206

0.10

0.5

0.7176 0.7828 0.8342 0.8713 0.8936 0.9011 0.8936 0.8713 0.8342 0.7828 0.7176 0.3619 0.3949 0.4212 0.4403 0.4519 0.4558 0.4519 0.4403 0.4212 0.3949 0.3619

0.6869 0.7556 0.8102 0.8498 0.8738 0.8818 0.8738 0.8498 0.8102 0.7556 0.6869 0.3283 0.3611 0.3873 0.4063 0.4179 0.4218 0.4179 0.4063 0.3873 0.3611 0.3283

0.20

1.00

0.6040 0.6591 0.7029 0.7348 0.7541 0.7606 0.7541 0.7348 0.7029 0.6591 0.6040 0.1542 0.1682 0.1794 0.1875 0.1925 0.1941 0.1925 0.1875 0.1794 0.1682 0.1542

0.5709 0.6280 0.6735 0.7067 0.7268 0.7336 0.7268 0.7067 0.6735 0.6280 0.5709 0.1305 0.1435 0.1540 0.1615 0.1661 0.1677 0.1661 0.1615 0.1540 0.1435 0.1305

表 8:取 ?x ? 0.1, ? ? 0.25, 用中心差商(二阶)代替微商计算

t ? 0.005,0.01,0.015,0.02,0.10,0.20,0.25,0.50,1.00 时刻的解析解与数值解
时间(t) 解析解 数值解 0.9275 0.9875 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9875 0.9275 0.8590 0.9296 0.9708 0.9902 时间(t) 解析解 数值解 0.8978 0.9648 0.9923 0.9992 1.0000 1.0000 1.0000 0.9992 0.9923 0.9648 0.8978 0.7175 0.7829 0.8345 0.8718 时间(t) 解析解 数值解 0.8764 0.9459 0.9818 0.9956 0.9993 0.9999 0.9993 0.9956 0.9818 0.9459 0.8764 0.6037 0.6589 0.7029 0.7348

0.005

0.9250 0.9841 0.9984 0.9999 1.0000 1.0000 1.0000 0.9999 0.9984 0.9841 0.9250 0.8585 0.9286 0.9695 0.9891

0.01

0.8965 0.9627 0.9905 0.9984 0.9998 1.0000 0.9998 0.9984 0.9905 0.9627 0.8965 0.7176 0.7828 0.8342 0.8713
- 14 -

0.015

0.8755 0.9444 0.9802 0.9945 0.9988 0.9996 0.9988 0.9945 0.9802 0.9444 0.8755 0.6040 0.6591 0.7029 0.7348

0.02

0.25

0.9967 0.9985 0.9967 0.9891 0.9695 0.9286 0.8585 0.5546 0.6052 0.6454 0.6747 0.6924 0.6984 0.6924 0.6747 0.6454 0.6052 0.5546

0.9974 0.9991 0.9974 0.9902 0.9708 0.9296 0.8590

0.10

0.5541 0.6048 0.6452 0.6745 0.6923 0.6983 0.6923 0.6745 0.6452 0.6048 0.5541

0.50

0.8936 0.9011 0.8936 0.8713 0.8342 0.7828 0.7176 0.3619 0.3949 0.4212 0.4403 0.4519 0.4558 0.4519 0.4403 0.4212 0.3949 0.3619

0.8942 0.9017 0.8942 0.8718 0.8345 0.7829 0.7175

0.20

0.3612 0.3942 0.4205 0.4396 0.4512 0.4551 0.4512 0.4396 0.4205 0.3942 0.3612

1.00

0.7541 0.7606 0.7541 0.7348 0.7029 0.6591 0.6040 0.1542 0.1682 0.1794 0.1875 0.1925 0.1941 0.1925 0.1875 0.1794 0.1682 0.1542

0.7542 0.7607 0.7542 0.7348 0.7029 0.6589 0.6037

0.1534 0.1674 0.1786 0.1867 0.1917 0.1933 0.1917 0.1867 0.1786 0.1674 0.1534

图 7:取 ?x ? 0.1, ? ? 0.25, 用向前(后)一阶差商代替微商计算

t ? 0.005,0.01,0.015,0.02,0.10,0.20,0.25,0.50,1.00 时刻的解析解与数值解

- 15 -

图 8:取 ?x ? 0.1, ? ? 0.25, 用中心差商(二阶)代替微商计算

t ? 0.005,0.01,0.015,0.02,0.10,0.20,0.25,0.50,1.00 时刻的解析解与数值解

例 4 数值求解边值问题 ? ut ? a ? x, t ? u x ? 0, x ? 0, t ? 0 ? ?1 当 0.2 ? x ? 0.4 ? ?u ( x, 0) ? ? 其它 ?0 ? ? u (0,t)=0, t?0 ? 其中 a( x, t ) ?
1 ? x2 , x ? 0, t ? 0 1 ? 2 xt ? 2 x 2 ? x 4

要求采用下列方法计算: (1) 迎风格式,分别取 ?x ? 0.01,0.02, t ? 0.1,0.5,1.0 时刻的数值解。 (2) Lax ? Friedrichs 格式完成(1)的计算。 (3) Lax ? Wendroff 格式完成(1)的计算。 注:该初边值问题的解析解为 t u ( x, t ) ? u ( x ? , 0) 1 ? x2 解:程序 function A=yingfenggeshi(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4); end end U=zeros(length(x),length(y));
- 16 -

for i=1:length(x) if x(i)>=0.2&x(i)<=0.4 U(i,1)=1; else U(i,1)=0; end end for j=1:length(y)-1 for i=1:length(x)-1 U(i+1,j+1)=(1-a(i+1,j+1)*la)*U(i+1,j)+a(i+1,j+1)*la*U(i,j); end end A=U(:,size(U,2)); function A=laxfriedrichs(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4); end end U=zeros(length(x),length(y)); for i=1:length(x) if x(i)>=0.2&x(i)<=0.4 U(i,1)=1; else U(i,1)=0; end end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=(1-a(i+1,j+1)*la)*U(i+2,j)/2+(1+a(i+1,j+1)*la)*U(i,j)/2; end end A=U(:,size(U,2)); function A=laxwendroff(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4); end end U=zeros(length(x),length(y)); for i=1:length(x) if x(i)>=0.2&x(i)<=0.4
- 17 -

U(i,1)=1; else U(i,1)=0; end end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)-a(i+1,j+1)*la*(U(i+2,j)-U(i,j))/2+a(i+1,j+1)^2*la^2*(U(i+2,j)-2 *U(i+1,j)+U(i,j))/2; end end A=U(:,size(U,2)); function u=jiexijie3(x,t) for i=1:length(x) if x(i)-t/(1+x(i)^2)>=0.2&x(i)-t/(1+x(i)^2)<=0.4 u(i)=1; else u(i)=0; end end u clc;%第 4 题第 1 问 clear; la=0.25;dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; x1=[0:dx1:1];x2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; A1=yingfenggeshi(x1,y1,la);u1=jiexijie3(x1,t1); A2=yingfenggeshi(x1,y2,la);u2=jiexijie3(x1,t2); A3=yingfenggeshi(x1,y3,la);u3=jiexijie3(x1,t3); A10=yingfenggeshi(x2,y10,la);u10=jiexijie3(x2,t1); A20=yingfenggeshi(x2,y20,la);u20=jiexijie3(x2,t2); A30=yingfenggeshi(x2,y30,la);u30=jiexijie3(x2,t3); subplot(121); line(x1,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x1,A2,'color','r','linestyle',':','linewidth',1.5); line(x1,A3,'color','r','linestyle',':','linewidth',1.5); line(x1,u1,'color','b','linewidth',1); line(x1,u2,'color','b','linewidth',1); line(x1,u3,'color','b','linewidth',1); title('例 4-1-yingfenggeshi');
- 18 -

subplot(122); line(x2,A10,'color','r','linestyle',':','linewidth',1.5); line(x2,A20,'color','r','linestyle',':','linewidth',1.5); line(x2,A30,'color','r','linestyle',':','linewidth',1.5); line(x2,u10,'color','b','linewidth',1); line(x2,u20,'color','b','linewidth',1); line(x2,u30,'color','b','linewidth',1); title('例 4-1-yingfenggeshi'); clc;%第 4 题第 2 问 clear; la=0.25;dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; X1=[0:dx1:1];X2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; x1=[0:dx1:(length(X1)+length(y1)-1)*dx1] x2=[0:dx1:(length(X1)+length(y2)-1)*dx1]; x3=[0:dx1:(length(X1)+length(y3)-1)*dx1]; x10=[0:dx2:(length(X2)+length(y10)-1)*dx2]; x20=[0:dx2:(length(X2)+length(y20)-1)*dx2]; x30=[0:dx2:(length(X2)+length(y30)-1)*dx2]; A1=laxfriedrichs(x1,y1,la);u1=jiexijie3(X1,t1); A2=laxfriedrichs(x2,y2,la);u2=jiexijie3(X1,t2); A3=laxfriedrichs(x3,y3,la);u3=jiexijie3(X1,t3); A10=laxfriedrichs(x10,y10,la);u10=jiexijie3(X2,t1); A20=laxfriedrichs(x20,y20,la);u20=jiexijie3(X2,t2); A30=laxfriedrichs(x30,y30,la);u30=jiexijie3(X2,t3); subplot(121); line(X1,A1(1:length(X1)),'color','r','linestyle',':','linewidth',1.5);hold on line(X1,A2(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,A3(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,u1(1:length(X1)),'color','b','linewidth',1); line(X1,u2(1:length(X1)),'color','b','linewidth',1); line(X1,u3(1:length(X1)),'color','b','linewidth',1);title('例 4-2-laxfriedrichs'); subplot(122); line(X2,A10(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A20(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A30(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,u10(1:length(X2)),'color','b','linewidth',1); line(X2,u20(1:length(X2)),'color','b','linewidth',1); line(X2,u30(1:length(X2)),'color','b','linewidth',1);title('例 4-2-laxfriedrichs'); clc;%第 4 题第 3 问 clear;la=0.25;
- 19 -

dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; X1=[0:dx1:1];X2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; x1=[0:dx1:(length(X1)+length(y1)-1)*dx1]; x2=[0:dx1:(length(X1)+length(y2)-1)*dx1]; x3=[0:dx1:(length(X1)+length(y3)-1)*dx1]; x10=[0:dx2:(length(X2)+length(y10)-1)*dx2]; x20=[0:dx2:(length(X2)+length(y20)-1)*dx2]; x30=[0:dx2:(length(X2)+length(y30)-1)*dx2]; A1=laxwendroff(x1,y1,la);u1=jiexijie3(X1,t1); A2=laxwendroff(x2,y2,la);u2=jiexijie3(X1,t2); A3=laxwendroff(x3,y3,la);u3=jiexijie3(X1,t3); A10=laxwendroff(x10,y10,la);u10=jiexijie3(X2,t1); A20=laxwendroff(x20,y20,la);u20=jiexijie3(X2,t2); A30=laxwendroff(x30,y30,la);u30=jiexijie3(X2,t3); subplot(121); line(X1,A1(1:length(X1)),'color','r','linestyle',':','linewidth',1.5);hold on line(X1,A2(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,A3(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,u1(1:length(X1)),'color','b','linewidth',1); line(X1,u2(1:length(X1)),'color','b','linewidth',1); line(X1,u3(1:length(X1)),'color','b','linewidth',1);title('例 4-3-laxwendroff'); subplot(122); line(X2,A10(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A20(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A30(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,u10(1:length(X2)),'color','b','linewidth',1); line(X2,u20(1:length(X2)),'color','b','linewidth',1); line(X2,u30(1:length(X2)),'color','b','linewidth',1);title('例 4-3-laxwendroff'); 运行结果:

- 20 -

图9:取 ?x ? 0.01, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与yingfenggeshi数值解

图10:取 ?x ? 0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与yingfenggeshi数值解

- 21 -

图11:取 ?x ? 0.01, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与laxfriedrichs数值解

图12:取 ?x ? 0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与laxfriedrichs数值解

- 22 -

图13:取 ?x ? 0.01, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与laxwendroff数值解

图 14:取 ?x ? 0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与 laxwendroff 数值解

例 5 数值求解边值问题
?ut ? a ? x, t ? u x ? 0, x ? 0, t ? 0 ? 2 2 ? u ( x, 0) ? exp(?10(4 x ? 1) ) ? u (0,t)=0, t?0 ?

1 ? x2 , x ? 0, t ? 0 其中 a( x, t ) ? 1 ? 2 xt ? 2 x 2 ? x 4

要求采用下列方法计算: (4) 迎风格式,分别取 ?x ? 0.01,0.02, t ? 0.1,0.5,1.0 时刻的数值解。
- 23 -

(5) Lax ? Friedrichs 格式完成(1)的计算。 (6) Lax ? Wendroff 格式完成(1)的计算。 注:该初边值问题的解析解为 t u ( x, t ) ? u ( x ? , 0) 1 ? x2 解:程序: function A=yingfenggeshi1(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4); end end U=zeros(length(x),length(y)); for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-1 U(i+1,j+1)=(1-a(i+1,j+1)*la)*U(i+1,j)+a(i+1,j+1)*la*U(i,j); end end A=U(:,size(U,2)); function A=laxfriedrichs1(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4); end end U=zeros(length(x),length(y)); for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=(1-a(i+1,j+1)*la)*U(i+2,j)/2+(1+a(i+1,j+1)*la)*U(i,j)/2; end end A=U(:,size(U,2)); function A=laxwendroff1(x,y,la) for i=1:length(x) for j=1:length(y) a(i,j)=(1+x(i)^2)/(1+2*x(i)*y(j)+2*x(i)^2+x(i)^4);
- 24 -

end end U=zeros(length(x),length(y)); for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)-a(i+1,j+1)*la*(U(i+2,j)-U(i,j))/2+a(i+1,j+1)^2*la^2*(U(i+2,j)-2 *U(i+1,j)+U(i,j))/2; end end A=U(:,size(U,2)); clc;%第 5 题第 1 问 clear; la=0.25;dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; x1=[0:dx1:1];x2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; A1=yingfenggeshi1(x1,y1,la);u1=jiexijie4(x1,t1); A2=yingfenggeshi1(x1,y2,la);u2=jiexijie4(x1,t2); A3=yingfenggeshi1(x1,y3,la);u3=jiexijie4(x1,t3); A10=yingfenggeshi1(x2,y10,la);u10=jiexijie4(x2,t1); A20=yingfenggeshi1(x2,y20,la);u20=jiexijie4(x2,t2); A30=yingfenggeshi1(x2,y30,la);u30=jiexijie4(x2,t3); subplot(121); line(x1,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x1,A2,'color','r','linestyle',':','linewidth',1.5); line(x1,A3,'color','r','linestyle',':','linewidth',1.5); line(x1,u1,'color','b','linewidth',1); line(x1,u2,'color','b','linewidth',1); line(x1,u3,'color','b','linewidth',1);title('例 5-1-yingfenggeshi'); subplot(122); line(x2,A10,'color','r','linestyle',':','linewidth',1.5); line(x2,A20,'color','r','linestyle',':','linewidth',1.5); line(x2,A30,'color','r','linestyle',':','linewidth',1.5); line(x2,u10,'color','b','linewidth',1); line(x2,u20,'color','b','linewidth',1); line(x2,u30,'color','b','linewidth',1);% title('例 5-1-yingfenggeshi'); clc;%第 5 题第 2 问 clear;
- 25 -

la=0.25;dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; X1=[0:dx1:1];X2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; x1=[0:dx1:(length(X1)+length(y1)-1)*dx1]; x2=[0:dx1:(length(X1)+length(y2)-1)*dx1]; x3=[0:dx1:(length(X1)+length(y3)-1)*dx1]; x10=[0:dx2:(length(X2)+length(y10)-1)*dx2]; x20=[0:dx2:(length(X2)+length(y20)-1)*dx2]; x30=[0:dx2:(length(X2)+length(y30)-1)*dx2]; A1=laxfriedrichs1(x1,y1,la);u1=jiexijie4(X1,t1); A2=laxfriedrichs1(x2,y2,la);u2=jiexijie4(X1,t2); A3=laxfriedrichs1(x3,y3,la);u3=jiexijie4(X1,t3); A10=laxfriedrichs1(x10,y10,la);u10=jiexijie4(X2,t1); A20=laxfriedrichs1(x20,y20,la);u20=jiexijie4(X2,t2); A30=laxfriedrichs1(x30,y30,la);u30=jiexijie4(X2,t3); subplot(121); line(X1,A1(1:length(X1)),'color','r','linestyle',':','linewidth',1.5);hold on line(X1,A2(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,A3(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,u1(1:length(X1)),'color','b','linewidth',1); line(X1,u2(1:length(X1)),'color','b','linewidth',1); line(X1,u3(1:length(X1)),'color','b','linewidth',1);title('例 5-2-laxfriedrichs'); subplot(122); line(X2,A10(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A20(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A30(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,u10(1:length(X2)),'color','b','linewidth',1); line(X2,u20(1:length(X2)),'color','b','linewidth',1); line(X2,u30(1:length(X2)),'color','b','linewidth',1);% title('例 5-2-laxfriedrichs'); clc;%第 5 题第 3 问 clear; la=0.25;dx1=0.01;dx2=0.02;t1=0.1;t2=0.5;t3=1.0; X1=[0:dx1:1];X2=[0:dx2:1]; y1=[0:dx1*la:t1];y2=[0:dx1*la:t2];y3=[0:dx1*la:t3]; y10=[0:dx2*la:t1];y20=[0:dx2*la:t2];y30=[0:dx2*la:t3]; x1=[0:dx1:(length(X1)+length(y1)-1)*dx1]; x2=[0:dx1:(length(X1)+length(y2)-1)*dx1]; x3=[0:dx1:(length(X1)+length(y3)-1)*dx1]; x10=[0:dx2:(length(X2)+length(y10)-1)*dx2]; x20=[0:dx2:(length(X2)+length(y20)-1)*dx2]; x30=[0:dx2:(length(X2)+length(y30)-1)*dx2]; A1=laxwendroff1(x1,y1,la);u1=jiexijie4(X1,t1);
- 26 -

A2=laxwendroff1(x2,y2,la);u2=jiexijie4(X1,t2); A3=laxwendroff1(x3,y3,la);u3=jiexijie4(X1,t3); A10=laxwendroff1(x10,y10,la);u10=jiexijie4(X2,t1); A20=laxwendroff1(x20,y20,la);u20=jiexijie4(X2,t2); A30=laxwendroff1(x30,y30,la);u30=jiexijie4(X2,t3); subplot(121); line(X1,A1(1:length(X1)),'color','r','linestyle',':','linewidth',1.5);hold on line(X1,A2(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,A3(1:length(X1)),'color','r','linestyle',':','linewidth',1.5); line(X1,u1(1:length(X1)),'color','b','linewidth',1); line(X1,u2(1:length(X1)),'color','b','linewidth',1); line(X1,u3(1:length(X1)),'color','b','linewidth',1); title('例 5-3-laxwendroff'); subplot(122); line(X2,A10(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A20(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,A30(1:length(X2)),'color','r','linestyle',':','linewidth',1.5); line(X2,u10(1:length(X2)),'color','b','linewidth',1); line(X2,u20(1:length(X2)),'color','b','linewidth',1); line(X2,u30(1:length(X2)),'color','b','linewidth',1);% title('例 5-3-laxwendroff');、 运行结果:

图 15:取 ?x ? 0.01,0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与 yingfenggesh 数值解

- 27 -

图 16:取 ?x ? 0.01,0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与 laxfriedrichs 数值解

图 17:取 ?x ? 0.01,0.02, ? ? 0.25, t ? 0.1,0.5,1.0, 时刻的解析解与 laxwendroff 数值解

1 2 ? ? ut ? ( u ) x ? 0, x ? R, t ? 0 例 6 问题 ? 2 2 2 ? ?u ( x, 0) ? exp(?10(4 x ? 1) ), x ? R

试用迎风格式, Lax ? Friedrichs 格式, Lax ? Wendroff 格式,计算该问题在
t ? 0.15, 0.3 时刻的数值解。

注:解析解为 u( x, t ) ? u0 ( x ? tu( x, t )) 解:程序:
- 28 -

function A=yingfenggeshi2(x,y,la) U=zeros(length(x),length(y)); for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-1 U(i+1,j+1)=U(i+1,j)-la*(U(i+1,j)^2/2-U(i,j)^2/2); end end A=U(:,size(U,2)); function A=laxfriedrichs2(x,y,la) U=zeros(length(x),length(y)); for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=(U(i,j)+U(i+2,j))/2-la*(U(i+2,j)^2/2-U(i,j)^2/2)/2; end end A=U(:,size(U,2)); function A=laxwendroff2(x,y,la) for i=1:length(x) U(i,1)=exp(-10*(4*x(i)^2-1)^2); end for j=1:length(y)-1 for i=1:length(x)-2 U(i+1,j+1)=U(i+1,j)-la*(U(i+2,j)^2-U(i,j)^2)/4+la^2*(U(i+1,j)+U(i+2,j))*(U(i+2,j)^2 -U(i+1,j)^2)/8-la^2*(U(i,j)+U(i+1,j))*(U(i+1,j)^2-U(i,j)^2)/8; end end A=U(:,size(U,2)); function u=jiexijie5(x,t) u(1)=0; for i=2:length(x) u(i)=exp(-10*(4*(x(i)-t*u(i-1))^2-1)^2); end clc;%第 6 题第 1 问
- 29 -

clear; la=0.3;t1=0.15;t2=0.3;x1=[0:0.01:1];x2=[0:0.02:1]; y1=[0:0.01*la:t1];y2=[0:0.01*la:t2]; y10=[0:0.02*la:t1];y20=[0:0.02*la:t2]; A1=yingfenggeshi2(x1,y1,la);u1=jiexijie5(x1,t1); A2=yingfenggeshi2(x1,y2,la);u2=jiexijie5(x1,t2); A10=yingfenggeshi2(x2,y10,la);u10=jiexijie5(x2,t1); A20=yingfenggeshi2(x2,y20,la);u20=jiexijie5(x2,t2); subplot(121) title('dx=0.01,t=0.15'); line(x1,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x1,u1,'color','b','linewidth',1); subplot(122) line(x1,A2,'color','r','linestyle',':','linewidth',1.5); line(x1,u2,'color','b','linewidth',1);title('dx=0.01,t=0.3'); figure; subplot(121) title('dx=0.02,t=0.15'); line(x2,A10,'color','r','linestyle',':','linewidth',1.5); line(x2,u10,'color','b','linewidth',1); subplot(122) line(x2,A20,'color','r','linestyle',':','linewidth',1.5); line(x2,u20,'color','b','linewidth',1);title('dx=0.02,t=0.3'); clc;%第 6 题第 2 问 clear;la=0.3;t1=0.15;t2=0.3; x1=[0:0.01:1];x2=[0:0.02:1]; y1=[0:0.01*la:t1];y2=[0:0.01*la:t2]; y10=[0:0.02*la:t1];y20=[0:0.02*la:t2]; A1=laxfriedrichs2(x1,y1,la);A2=laxfriedrichs2(x1,y2,la); A10=laxfriedrichs2(x2,y10,la);A20=laxfriedrichs2(x2,y20,la); u1=jiexijie5(x1,t1);u2=jiexijie5(x1,t2); u10=jiexijie5(x2,t1);u20=jiexijie5(x2,t2); subplot(121) title('dx=0.01,t=0.15'); line(x1,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x1,u1,'color','b','linewidth',1); subplot(122) line(x1,A2,'color','r','linestyle',':','linewidth',1.5); line(x1,u2,'color','b','linewidth',1); title('dx=0.01,t=0.3'); figure; subplot(121) title('dx=0.02,t=0.15');
- 30 -

line(x2,A10,'color','r','linestyle',':','linewidth',1.5); line(x2,u10,'color','b','linewidth',1); subplot(122) line(x2,A20,'color','r','linestyle',':','linewidth',1.5); line(x2,u20,'color','b','linewidth',1); title('dx=0.02,t=0.3'); clc;%第 6 题第 3 问 clear; la=0.3;t1=0.15;t2=0.3; x1=[0:0.01:1];x2=[0:0.02:1]; y1=[0:0.01*la:t1];y2=[0:0.01*la:t2]; y10=[0:0.02*la:t1];y20=[0:0.02*la:t2]; A1=laxwendroff2(x1,y1,la);A2=laxwendroff2(x1,y2,la); A10=laxwendroff2(x2,y10,la);A20=laxwendroff2(x2,y20,la); u1=jiexijie5(x1,t1);u2=jiexijie5(x1,t2); u10=jiexijie5(x2,t1);u20=jiexijie5(x2,t2); subplot(121) title('dx=0.01,t=0.15'); line(x1,A1,'color','r','linestyle',':','linewidth',1.5);hold on line(x1,u1,'color','b','linewidth',1); subplot(122) line(x1,A2,'color','r','linestyle',':','linewidth',1.5); line(x1,u2,'color','b','linewidth',1); title('dx=0.01,t=0.3'); figure; subplot(121) title('dx=0.02,t=0.15'); line(x2,A10,'color','r','linestyle',':','linewidth',1.5); line(x2,u10,'color','b','linewidth',1); subplot(122) line(x2,A20,'color','r','linestyle',':','linewidth',1.5); line(x2,u20,'color','b','linewidth',1); title('dx=0.02,t=0.3'); 运行结果:

- 31 -

图 18:取 ?x ? 0.01, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 yingfenggesh 数值解

图 19:取 ?x ? 0.02, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 yingfenggesh 数值解

- 32 -

图 20:取 ?x ? 0.01, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 laxfriedrichs 数值解

图 21:取 ?x ? 0.02, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 laxfriedrichs 数值解

- 33 -

图 22:取 ?x ? 0.01, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 laxwendroff 数值解

图 23:取 ?x ? 0.02, ? ? 0.3 t ? 0.15, 0.3 时刻的解析解与 laxwendroff 数值解

全部程序完成

- 34 -


相关文章:
偏微分方程数值解法考试题
偏微分方程数值解法考试题_理学_高等教育_教育专区。2014—2015 学年第二学期...的常用的差分方法进行归纳总结,分析它们各自的优缺点,算出其精确解,用前面的...
偏微分方程数值解题目
偏微分方程数值解题目_数学_自然科学_专业资料。题目一、用有限元法求下列边值问题的数值解: ? ?2 ?x y = 2sin( ),0 < x < 1, ?-y ''+ 4 2 ?...
偏微分方程数值解试题参考答案
偏微分方程数值解试题参考答案 - 偏微分方程数值解 1 ( Ax , x) ? (b, x) ( x ? R n ) ,证明下 2 一(10 分) 、设矩阵 A 对称正定,定义 J (...
偏微分方程数值解习题解答案
偏微分方程数值解习题解答案_理学_高等教育_教育专区。 文档贡献者 滕卜卜 ...偏微分方程数值解 13页 1下载券 偏微分方程数值解(试题) 7页 2下载券 喜欢...
2014-2015学年偏微分方程数值解法试题
2014-2015学年偏微分方程数值解法试题 - 编号 浙江理工大学考试命题稿( A (2014 /2015 学年 第 1 开课学院: 机械与自动控制学院 卷) 学期) 开课年级 课程...
偏微分方程数值解习题解答案
偏微分方程数值解习题解答案 - q1 q2 2 第二章 第三章 第四章 第五章 第六章 3页 A 1 Q 2 13 页 q3 a3 q4 a4 q5 2...
偏微分方程数值习题解答
偏微分方程数值习题解答 - 李微分方程数值解习题解答 1-1 如果? ' (0) = 0 ,则称 x0 是 J (x) 的 对称(不必正定) 驻点(或稳定点) , 驻点(或稳定...
偏微分方程数值解例题答案
偏微分方程数值解例题答案 - 二、改进的 Euler 方法 梯形方法的迭代公式(1.10)比 Euler 方法精度高,但其计算较复杂,在应用公式(1.10)进行 计算时,每迭代一次,...
偏微分方程数值解期末试题及答案
偏微分方程数值解试题(06B) 偏微分方程数值解试题(06B) 参考答案与评分标准信息与计算科学专业 一( 10 分)、设矩阵 A 对称,定义 J ( x) = 1 ( Ax, x)...
偏微分方程数值解复习题(2013硕士)
偏微分方程数值解复习题(2013硕士) - 偏微分方程数值解期末复习(2012 硕士) 一、 考题类型 本次试卷共六道题目, 题型及其所占比例分别为: 填空题 20%; 计算...
更多相关标签: