微分方程式を行列で差分化して解く その２
前回からの続きです。 前回は以下
今回は計算にJulia言語を使用しています。この文章自体も一部Weave.jlを使って書いています。
まず使用するパッケージを読み込んでおきます。
using Gadfly
using LinearAlgebra
前回に引き続き以下の微分方程式を行列を使用して解きます。
前回は以下のようになりました。
A = [1 0 0 0;-1 1 0 0;0 -1 1 0;0 0 -1 1]
b = [0 3 12 27]'
y = A\b
println(y)
[0.0; 3.0; 15.0; 42.0]
微分方程式の解析解は
ですが今回のyはこのグラフ上に乗るでしょうか？一応プロットしてみます。
x1 = 0:0.1:3
x2 = 0:1:3
y_real = x1.^3
plot(
layer(x = x1, y = y_real,Geom.line,Theme(default_color=colorant"red")),
layer(x = x2, y = y)
)
………全然違いますね。これは前進差分の⊿x = 1なので大きく取りすぎたものと考えられます。 もうすこし差分を小さくして、例えば⊿x = 0.01 として計算してみたいと思います。
その場合、数式としては以下のようになりますね。
これを解いてみます。
dv = ones(301)
ev = -1 *ones(300)
A_new = Bidiagonal(dv,ev, :L)
k = 0:0.01:3
b_new = reshape(0.01*3*k.^2,301,1)
y_new = A_new \ b_new
println(y_new)
[0.0; 3.0e-6; 1.5e-5; 4.2e-5; 9.0e-5; 0.000165; 0.000273; 0.00042; 0.000612
; 0.000855; 0.001155; 0.001518; 0.00195; 0.002457; 0.003045; 0.00372; 0.004
488; 0.005355; 0.006327; 0.00741; 0.00861; 0.009933; 0.011385; 0.012972; 0.
0147; 0.016575; 0.018603; 0.02079; 0.023142; 0.025665; 0.028365; 0.031248;
0.03432; 0.037587; 0.041055; 0.04473; 0.048618; 0.052725; 0.057057; 0.06162
; 0.06642; 0.071463; 0.076755; 0.082302; 0.08811; 0.094185; 0.100533; 0.107
16; 0.114072; 0.121275; 0.128775; 0.136578; 0.14469; 0.153117; 0.161865; 0.
17094; 0.180348; 0.190095; 0.200187; 0.21063; 0.22143; 0.232593; 0.244125;
0.256032; 0.26832; 0.280995; 0.294063; 0.30753; 0.321402; 0.335685; 0.35038
5; 0.365508; 0.38106; 0.397047; 0.413475; 0.43035; 0.447678; 0.465465; 0.48
3717; 0.50244; 0.52164; 0.541323; 0.561495; 0.582162; 0.60333; 0.625005; 0.
647193; 0.6699; 0.693132; 0.716895; 0.741195; 0.766038; 0.79143; 0.817377;
0.843885; 0.87096; 0.898608; 0.926835; 0.955647; 0.98505; 1.01505; 1.04565;
1.07686; 1.10869; 1.14114; 1.17421; 1.20792; 1.24227; 1.27726; 1.3129; 1.3
492; 1.38617; 1.4238; 1.46211; 1.5011; 1.54077; 1.58114; 1.6222; 1.66398; 1
.70646; 1.74966; 1.79358; 1.83823; 1.88362; 1.92975; 1.97662; 2.02425; 2.07
264; 2.12179; 2.17171; 2.22241; 2.2739; 2.32617; 2.37924; 2.4331; 2.48778;
2.54327; 2.59957; 2.65671; 2.71467; 2.77347; 2.83311; 2.8936; 2.95495; 3.01
716; 3.08024; 3.14418; 3.20901; 3.27472; 3.34133; 3.40883; 3.47723; 3.54654
; 3.61677; 3.68792; 3.75999; 3.833; 3.90695; 3.98184; 4.05768; 4.13448; 4.2
1224; 4.29098; 4.37068; 4.45137; 4.53305; 4.61571; 4.69938; 4.78405; 4.8697
4; 4.95644; 5.04416; 5.13291; 5.2227; 5.31353; 5.4054; 5.49833; 5.59232; 5.
68737; 5.78349; 5.88069; 5.97897; 6.07835; 6.17881; 6.28038; 6.38306; 6.486
84; 6.59175; 6.69778; 6.80495; 6.91325; 7.02269; 7.13328; 7.24503; 7.35794;
7.47201; 7.58726; 7.70369; 7.8213; 7.9401; 8.0601; 8.1813; 8.30372; 8.4273
4; 8.55219; 8.67827; 8.80557; 8.93412; 9.06391; 9.19496; 9.32726; 9.46082;
9.59565; 9.73176; 9.86915; 10.0078; 10.1478; 10.2891; 10.4316; 10.5755; 10.
7207; 10.8672; 11.0151; 11.1643; 11.3148; 11.4667; 11.6199; 11.7745; 11.930
4; 12.0878; 12.2465; 12.4065; 12.568; 12.7309; 12.8952; 13.0608; 13.2279; 1
3.3964; 13.5664; 13.7377; 13.9105; 14.0848; 14.2605; 14.4376; 14.6162; 14.7
963; 14.9778; 15.1609; 15.3454; 15.5314; 15.7189; 15.9079; 16.0984; 16.2904
; 16.484; 16.679; 16.8756; 17.0738; 17.2735; 17.4747; 17.6775; 17.8819; 18.
0878; 18.2953; 18.5044; 18.7151; 18.9274; 19.1412; 19.3567; 19.5738; 19.792
5; 20.0128; 20.2348; 20.4583; 20.6836; 20.9104; 21.139; 21.3692; 21.601; 21
.8345; 22.0697; 22.3066; 22.5452; 22.7855; 23.0274; 23.2711; 23.5165; 23.76
36; 24.0124; 24.263; 24.5153; 24.7693; 25.0251; 25.2827; 25.542; 25.8031; 2
6.0659; 26.3305; 26.5969; 26.8651; 27.1351]
これを再びグラフにしてみます。
x = 0:0.01:3
y_real = x.^3
plot(
layer(x = x, y = y_real,Geom.line,Theme(default_color=colorant"red")),
layer(x = x, y = y_new)
)
今回は解析解とよく一致していますね。 ⊿xを細かく取ることがやはりポイントでした。
