科学计算课程设计-New

采用欧拉法
dy f ( x, y ) dx y ( x0 ) y0
y0 y ( x0 ) yi 1 yi hf ( xi , yi ) i 0...n 1
yi 1 yi hf ( xi , yi )
初始条件： t
0
时：I 0
n
dQ dt
n 1 t
1 dI I n 1 I n 2 dt
n
dI dt
n 1 t
由欧拉法求
1. 先由欧拉法计算
Q 'n1 Qn I n t , I 'n 1 I n ( Va Qn R I n )t L LC L
2.M-C方法计算定积分
g ( x) sin x / x 是[0, 1]上的连续函数,
求积分：
Ⅰ g ( x)dx
0
1
解：积分Ⅰ等于求图中曲线 g ( x) 下面曲线的面积。 向图中所示单位正方形内均匀投点(ξ,η),则该随机点落 入曲线g(x)下面的概率：
y 1 y=g(x) G 0 1 x
VR VL VC Va
由于： VR IR
dI VL L dt
Q VC C
dI Q IR L Va dt C
dQ dt I dI Va Q R I dt L LC L
d 2I dI I dV L 2 R a dt dt C dt
t 1
Va Qn R Qn1 Qn I n t , I n 1 I n ( I n )t L LC L
取
Q1 Q0 I 0 t 1 0 1 1 Va Q0 R Va Va 1 1 I1 I 0 ( I 0 )t 0 ( ) 1 L LC L L LC L LC
Pr y g ( x) g ( x)dx I
1 0
步骤： （1）产生两个随机数 i ， i。 （2）检验（ x i， i ）是否落在G内？是，则满足， 记下m=m+1。 设在M次试验中，共有m个随机点满足上述 （3） 不等式， 则 m Ⅰ M 即为积分近似值。
y
x y

3. M-C方法产生给定分布的随机变量
一、欧拉法/改进欧拉法-RCL回路问题
RLC串联电路: 电路由电阻 R、电感 L、电 容 C、开关 S 和外加电压Va, 组成。 求：当 开关S合上以后电路中的电流将随时间发生 变化 ？ 令R=L=C=1.0

பைடு நூலகம்
Va 5sin( t / 2)
分析：根据基尔霍夫第二定律，电路中的外电压等于电路内 各元件的电压之和
例：求具有分布密度
f ( x) 2 x
的随机抽样。
0 x 1
解： 设 Ri 是[0，1]区间上均匀分布的随机数，则
R 2 xdx 2 xdx x
0
x
x
2
xi Ri
x i 即为具有[0, 1]上且满足
f (x) 分布的随机数。
设计思路：
0 n (总次数); 0 N(次数); 0 M(命中数)
n n n (改变抽样总次数)
N=N+1
1 rand(1.0); 2 rand(1.0)
no (不计算)
1 2 1 2 1 (1 ) ( 2 ) ? 2 2 4
yes M=M+1 (命中数)
梯形公式
yi(0) yi hf ( xi , yi ) 1 ( k 1) h yi 1 yi ( f ( xi , yi ) f ( xi 1 , yi(k1) )) 2
改进梯形公式
tn 1 tn t
1 dQ Qn 1 Qn 2 dt
N n?
yes
no
n n n
4M N
no
0 ?
yes Print
其中
n1 n
小结：从此例可以看出用M-C法计算所需步骤： 1．把一个物理问题转化为一个数学问题; 2．这个数学问题一定包括随机数 (或随机抽样值); 3．产生随机数，并作试验。

利用[0,1]区间上的均匀分布随机数可以产生具有给定分 布的随机变量数列。
若随机变量具有分布密度f(x)，则随机变量F(x)的分布 就是区间[0, 1]上的均匀分布.
F ( x) f ( x)dx

x
因此，若Ri是[0, 1]中均匀分布的随机数，则方程
Ri f ( x)dx

2. 按照改进欧拉法计算
Qn 1 Qn 1 I n I n' 1 t 2
Va Q 'n 1 R 1 Va Qn R I n 1 I n ( In ) ( I 'n 1 ) t 2 L LC L L LC L
源程序： program main implicit none integer n, t, dt real r,l,c,q1(0:100),q2(0:100),i1(0:100),i2(0:100),qq(0:100),ii(0:100),va,pi parameter (pi=3.14159) !定义一个常量－圆周率 t=0 dt=1 r=1.0 l=1.0 c=1.0 q1(0)=1 q2(0)=1 i1(0)=0 i2(0)=0 open(unit=10,file="14.txt") write(10,"(a5,5x,a5,10x,a5)")"t","q2","i2" write(10,"(i5,5x,f10.5,5x,f10.5)")t,q2(0),i2(0) do n=0,99,1 va=5*sin(0.5*pi*t) qq(n)=i2(n) !qq(n)代表对电荷微分，ii(n)代表对电流的微分 ii(n)=(va-i2(n)*r-q2(n)/c)/l q1(n+1)=q2(n)+qq(n)*dt i1(n+1)=i2(n)+ii(n)*dt qq(n+1)=i1(n+1) ii(n+1)=(va-i1(n+1)*r-q1(n+1)/c)/l q2(n+1)=q2(n)+0.5*(qq(n)+qq(n+1))*dt !注意该式中的0.5不要写成1/2 i2(n+1)=i2(n)+0.5*(ii(n)+ii(n+1))*dt t=t+dt write(10,"(i5,5x,f10.5,5x,f10.5)")t,q2(n+1),i2(n+1) end do stop end
0
Q0 1
dQ dt I dI Va Q R I dt L LC L
tn 1 tn t dQ Qn 1 Qn t Qn I n t dt n
Va Qn R dI I n 1 I n t I n ( I n )t dt n L LC L
科学计算课程设计
武汉大学 机电学院
课程目的
巩固《计算方法》课程理论知识 学会利用相关计算方法求解实际物理问题 学习Monte Carlo算法 算例计算

课程要求
欧拉法/改进欧拉法计算RCL回路 Monte Carlo方法计算圆周率/定积分/随机 抽样

完成实验报告，并上交，作为最终成绩。
xi
的解xi 就是所求的具有分布密度为f(x)的随机抽样。
Ri f ( x)dx

xi
故，只要 F (x) 的反函数 F 1 ( x) 存在，
由[0, 1]区间均匀分布的随机数 。
xi F ( Ri ).
即解方程
1
Ri ,
求解
F ( x) R 就可得到分布函数为 f (x) 的随机抽样 xi 。
i
4
2
i(A)
0
-2
-4
0
20
40
60
80
100
t(s)
二、Monte Carlo方法

1. 在单位正方形内，有一内切园。若将针 均匀地投入正方形内，则命中圆内的概率 为：
P

4 1


4
y
如果投针N次，有M次命中圆内，
M 1 2 P 面积 ( ) N 2 4

x
4M 4P N
Q2 Q1 I1t Va Q0 R I 2 I1 ( I1 )t L LC L
……

改进欧拉法
yi(0) yi hf ( xi , yi ) 1 ( k 1) k yi 1 yi hf ( xi 1 , yi(1) )
y0 y ( x0 ) h yi 1 yi 2 ( f ( xi , yi ) f ( xi 1 , yi 1 ))