顯示具有 Julia 數值分析 標籤的文章。顯示所有文章

2019年3月7日星期四

Julia語言例題1-1 若已知下面二點請問 x=1.5時則f(1.5)= ??

Julia語言例題1-1 若已知下面二點

x y=f(x)
=================
1.0 0.000
2.0 0.693
請問 x=1.5時則f(1.5)= ??

程式
using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
#
# implements the interpolation algorithm of Newton
#
# ON ENTRY :
# x abscisses, given as a column vector;
# f ordinates, given as a column vector;
# xa point where to evaluate the interpolating
# polynomial through (x[i],f[i]).
#
# ON RETURN :
# d divided differences, computed from and f;
# p value of the interpolating polynomial at xa.
#
# EXAMPLE :
n = length(x)
tmp2=0.0
for k=1:n
tmp1=1.0
for i=1:n
if (i != k)
tmp1 *= (xa-x[i]) / (x[k]-x[i])
end
end
tmp2=tmp2+tmp1*f[k]
end

return tmp2
end

x = [1.0 , 2.0 ]
f = [0.0 , 0.693 ]
xa = 1.5
result1 = lagrange(x,f,xa)
println("Lagrange 內插法理則 ")
println("x= " , x)
println("f(x)= " , f)

s = @sprintf("Pn(x)=%0.5f" , result1 )
println(s)
s = @sprintf("f(xa)=%0.5f" , log(xa) )
println(s)
s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
println(s)

輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0]
f(x)= [0.0, 0.693]
Pn(x)=0.34650
f(xa)=0.40547
誤差 =0.05897

2019年3月5日星期二

Julia語言例題1-5 已知函數f(x)=ln(x)請使用Lagrange 內插法求f(1.5) ,f(2.5) f(3.5)的P(x)與f(x) 分別用x=1.5 ,2.5 ,3.5 計算f(x)的e(x)誤差值

Julia語言例題1-5 已知函數f(x)=ln(x)

xi     f(xi)

==============

1.0    0.000

2.0    0.693

3.0    1.099

4.0    1.386

請使用Lagrange 內插法 求f(1.5) ,f(2.5) f(3.5)的P(x)與f(x)

分別用x=1.5 ,2.5 ,3.5 計算f(x)的e(x)誤差值

using Printf

function lagrange(x::Array{Float64,1},f::Array{Float64,1},xa::Float64)
    #
    # implements the interpolation algorithm of Newton
    #
    # ON ENTRY :
    # x abscisses, given as a column vector;
    # f ordinates, given as a column vector;
    # xa point where to evaluate the interpolating
    # polynomial through (x[i],f[i]).
    #
    # ON RETURN :
    # d divided differences, computed from and f;
    # p value of the interpolating polynomial at xa.
    #
    # EXAMPLE :
    n = length(x)
    tmp2=0.0
    for k=1:n
        tmp1=1.0
        for i=1:n
            if (i != k)
                tmp1 *= (xa-x[i]) / (x[k]-x[i])
            end    
        end
        tmp2=tmp2+tmp1*f[k]    
    end    
        
    return tmp2
end


function er(xa::Float64,x0::Float64,x1::Float64,x2::Float64,x3::Float64)
    tmp1=0.0
    tmp2=0.0
    sum=1.0
    xm=(x[1]+x[length(x)])/2
    #println(xm)
    for i=1:length(x)
        sum=sum*i
    end 
    #println(sum)
    tmp1=(xa-x0)*(xa-x1)*(xa-x2)*(xa-x3)/sum
    #println(tmp1)
    tmp2=tmp1*(-6/xm^4)
    return tmp2
end

x= [1.0 ,2.0 , 3.0 , 4.0] 
f= [0.0,0.693,1.099,1.386]
xb = [ 1.5 , 2.5 ,3.5]

for i = 1:length(xb)
    xa=xb[i]
    result1 = lagrange(x,f,xa)
    println("Lagrange 內插法理則 ")
    println("x=    " , x)
    println("f(x)= " , f)
    println("xa=   " , xa)
    
   
    s = @sprintf("Pn(x)=%0.5f" , result1 )
    println(s)
    s = @sprintf("f(xa)=%0.5f" , log(xa) )
    println(s)
    s = @sprintf("誤差 =%0.5f" , abs(log(xa)-result1) )
    println(s)
    println()

end

for i=1:length(xb)
    result2=0.0
    xa=xb[i]
    x0=x[1]
    x1=x[2]
    x2=x[3]
    x3=x[4]
    result2= er(xa, x0 ,x1 , x2 , x3)
    s=@sprintf("x= %0.2f 誤差e(x)=%0.5f",xa ,result2)
    println(s)
end

輸出畫面

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   1.5
Pn(x)=0.39287
f(xa)=0.40547
誤差 =0.01259

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   2.5
Pn(x)=0.92138
f(xa)=0.91629
誤差 =0.00508

Lagrange 內插法理則 
x=    [1.0, 2.0, 3.0, 4.0]
f(x)= [0.0, 0.693, 1.099, 1.386]
xa=   3.5
Pn(x)=1.24687
f(xa)=1.25276
誤差 =0.00589

x= 1.50 誤差e(x)=0.00600
x= 2.50 誤差e(x)=-0.00360
x= 3.50 誤差e(x)=0.00600

2019年3月3日星期日

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 n-4 並比較誤差值

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 n=4 並比較誤差值

#==================================================
/* ex4-7-B.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= x exp(y)
return (x*exp(y))
end

function C(x::Float64) #// c(x)= 0.0
return (0.0)
end

function D(x::Float64) #// d(x)= x
return (x)
end

function gy(n::Int64)
sum1=0.0
sum2=0.0
for i=0:n
for j=1:n-1
if(i%2==0)
sum2=sum2+F(x[i+1],y[i+1][j+1])
else
sum1=sum1+F(x[i+1],y[i+1][j+1])
end
end
g1[i+1]=(hy[i+1]/3)*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum2+4*sum1)
println(i+1,"----",g1[i+1])
sum1=0.0
sum2=0.0
end
return g1
end

s=@sprintf("辛普森積分計算雙重積分 n=2 ")
println(s)

x= [0.0 for i=1:5 ]
g1=[0.0 for i=1:5 ]
g2=[0.0 for i=1:5 ]
hy=[0.0 for i=1:5 ]

f(i) = [0.0 for i=1:5]
y= f.([0.0 for i=1:5])

sum1=0.0
sum2=0.0

n=2
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
if(i%2==0)
sum2=sum2+g2[i+1]
else
sum1=sum1+g2[i+1]
end

end

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

s=@sprintf("梯形積分計算雙重積分 n=4 ")
println(s)

x= [0.0 for i=1:10 ]
g1=[0.0 for i=1:10 ]
g2=[0.0 for i=1:10 ]
hy=[0.0 for i=1:10 ]

f(i) = [0.0 for i=1:10]
y= f.([0.0 for i=1:10])

sum1=0.0
sum2=0.0

n=4
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

#println(x)
#println(hy)
#println(y)

g2=gy(n)
sum1=0.0
sum2=0.0

println("\n\n")
for i=1:n-1
if(i%2==0)
sum2=sum2+g2[i+1]
else
sum1=sum1+g2[i+1]
end

end

ts= (hx/3) * (g1[1] + g1[n+1] + 4*sum1 + 2*sum2)
s=@sprintf("辛普森積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

辛普森積分計算雙重積分 n=2 
1----0.0
2----0.3243676223937956
3----1.16928739497655



辛普森積分計算雙重積分結果 T2=0.411126

實際值=0.500000

誤差值=0.088874

梯形積分計算雙重積分 n=4 
1----0.0
2----0.0828099899078667
3----0.21652191332178472
4----0.9741611859661217
5----1.151481269705011



辛普森積分計算雙重積分結果 T4=0.484367

實際值=0.500000

誤差值=0.015633

2019年3月2日星期六

Julia語言例題EX4-7-A利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 與 n=4 二者並比較誤差值

Julia語言例題EX4-7-A利用梯形法求雙重積分
f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x
取 n=2 與 n=4 二者並比較誤差值

#==================================================
/* ex4-7-A.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= x exp(y)
return (x*exp(y))
end

function C(x::Float64) #// c(x)= 0.0
return (0.0)
end

function D(x::Float64) #// d(x)= x
return (x)
end

function gy(n::Int64)
sum=0.0;
for i=0:n
for j=1:n-1
sum=sum+F(x[i+1],y[i+1][j+1])
#println(y[i][j])
end
g1[i+1]=(0.5*hy[i+1])*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
println(i+1,"----",g1[i+1])
sum=0.0;
end
return g1
end

s=@sprintf("梯形積分計算雙重積分 n=2 ")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=2
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
sum1=sum1+g2[i+1];
#println(g2[i+1])
end

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

s=@sprintf("梯形積分計算雙重積分 n=4 ")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=4
a=0.0
b=1.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
sum1=sum1+g2[i+1];
#println(g2[i+1])
end

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(1/2))
println(s)
tn=abs( 1/2 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

梯形積分計算雙重積分 n=2 
1----0.0
2----0.3260482565047257
3----1.7539310924648253



梯形積分計算雙重積分結果 T2=0.601507

實際值=0.500000

誤差值=0.101507

梯形積分計算雙重積分 n=4 
1----0.0
2----0.07102946671483652
3----0.3247828699826771
4----0.8402029213086307
5----1.7272219045575166



梯形積分計算雙重積分結果 T4=0.524907

實際值=0.500000

誤差值=0.024907

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

#==================================================
/* ex4-6-C.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= x^2 * y
return (x*x*y)
end

function C(x::Float64) #// c(x)= 0.0
return (0.0)
end

function D(x::Float64) #// d(x)= 1.0
return (1.0)
end

function gy(n::Int64)
sum=0.0;
for i=0:n
for j=1:n-1
sum=sum+F(x[i+1],y[i+1][j+1])
#println(y[i][j])
end
g1[i+1]=(0.5*hy[i+1])*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
println(i+1,"----",g1[i+1])
sum=0.0;
end
return g1
end

s=@sprintf("梯形積分計算雙重積分")
println(s)

x= [0.0 for i=1:20 ]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=10
a=1.0
b=2.0

hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

g2=gy(n)
sum1=0.0
println("\n\n")
for i=1:n-1
sum1=sum1+g2[i+1];
#println(g2[i+1])
end

ts= (hx/2) * (g1[1] + g1[n+1] + 2*sum1)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%0.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%0.6lf\n",(7/6))
println(s)
tn=abs( 7/6 - ts )
s=@sprintf("誤差值=%0.6lf\n",tn )
println(s)

輸出畫面

梯形積分計算雙重積分
1----0.5000000000000001
2----0.6050000000000002
3----0.7200000000000001
4----0.8450000000000002
5----0.98
6----1.125
7----1.2800000000000005
8----1.4450000000000005
9----1.6200000000000003
10----1.8050000000000002
11----2.0000000000000004



梯形積分計算雙重積分結果 T10=1.167500

實際值=1.166667

誤差值=0.000833

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

Julia語言例題4-6 雙重定積分

b d

\ell =\iint _{T}f(x,y)\,dx\,dy

a c
a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

#==================================================
/* ex4-6.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= 4xy
return (x*x*y)
end

x=[0.0 for i=1:20]
g1=[0.0 for i=1:20 ]
y=[0.0 for i=1:20 ]

sum=0.0
n=2
a=1.0
b=2.0
hx=(b-a)/n
ax=a
for i=1:n+1
x[i]=ax
ax=ax+hx
end

c=0.0
d=1.0
hy=(d-c)/n
ay=c
for i=1:n+1
y[i]=ay
ay+=hy
end

for i=1:n+1
s=@sprintf("i=%2d , x[%1d]=%0.3f , y[%1d]=%0.3f ",i,i,x[i],i,y[i])
println(s)
end

println("\n梯形理則 (trapezoidal Rule) ")
Gx0=0.5*hy*(F(x[1],y[1]) + 2*F(x[1],y[2]) +F(x[1],y[3]))
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=0.5*hy*(F(x[2],y[1])+2* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=0.5*hy*(F(x[3],y[1])+2* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)

Tn=0.5*hx*(Gx0+2*Gx1+Gx2)
s=@sprintf("0.5*hy*(Gx0+2*Gx1+Gx2)= %0.5f",Tn)
println(s)

#===============================================#
println("\n辛普森理則（Simpson's Rule)")
Gx0=(1/3)*hy*( F(x[1],y[1])+4* F(x[1],y[2])+ F(x[1],y[3]) )
s=@sprintf("Gx0=%0.5f",Gx0)
println(s)

Gx1=(1/3)*hy*( F(x[2],y[1])+4* F(x[2],y[2])+ F(x[2],y[3]) )
s=@sprintf("Gx1=%0.5f",Gx1)
println(s)

Gx2=(1/3)*hy*( F(x[3],y[1])+4* F(x[3],y[2])+ F(x[3],y[3]) )
s=@sprintf("Gx2=%0.5f",Gx2)
println(s)
Tn=(1/3)*hy*(Gx0+4*Gx1+Gx2)
s=@sprintf("(1/3)*hy*(Gx0+4*Gx1+Gx2)= %0.5f", Tn)
println(s)
#===============================================#

輸出畫面

i= 1 , x[1]=1.000 , y[1]=0.000 
i= 2 , x[2]=1.500 , y[2]=0.500 
i= 3 , x[3]=2.000 , y[3]=1.000 

梯形理則 (trapezoidal Rule) 
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
0.5*hy*(Gx0+2*Gx1+Gx2)= 1.18750

辛普森理則（Simpson's Rule)
Gx0=0.50000
Gx1=1.12500
Gx2=2.00000
(1/3)*hy*(Gx0+4*Gx1+Gx2)= 1.16667

2019年3月1日星期五

Julia語言範例EX4-8 梯行積分法雙重積分計算

Julia語言範例EX4-8 梯行積分法雙重積分計算
n=10

b d(x)

\ell =\iint _{T}f(x,y)\,dx\,dy

a c(x)

a=0 , b=1 , c(x)=0 , d(x)=1-x^2 f(x,y)=4xy

#==================================================
/* ex4-8.jl based on Trapezoidal Rule to
* compute the double integral.
*/
===================================================#
using Printf

function F(x::Float64, y::Float64) #// f(x,y)= 4xy
return (4*x*y)
end

function C(x::Float64) #// c(x)=0
return (0.0)
end
function D(x::Float64) #// d(x)= 1-x^2
return (1-x*x)
end

function gy(n::Int64)
sum=0.0;
for i=0:n
for j=1:n-1
sum=sum+F(x[i+1],y[i+1][j+1])
#println(y[i][j])
end
g1[i+1]=(hy[i+1]/2.0)*(F(x[i+1],y[i+1][1])+F(x[i+1],y[i+1][n+1])+2*sum)
sum=0.0;
end
return g1
end

x=[0.0 for i=1:20]
g1=[0.0 for i=1:20 ]
g2=[0.0 for i=1:20 ]
hy=[0.0 for i=1:20 ]

f(i) = [0.0 for i=1:20]
y= f.([0.0 for i=1:20])

sum=0.0
n=10
a=0.0
b=1.0
hx=(b-a)/n
ts=0.0
for i=0:n
x[i+1]=a+i*hx
hy[i+1]=(D(x[i+1])-C(x[i+1]))/n
for j=0:n
y[i+1][j+1]=C(x[i+1])+j*hy[i+1]
#println(y)
end
end

g2=gy(n)
println("\n\n")
for i=1:n-1
sum=sum+g1[i+1];
end

ts=(hx/2.0)*(g1[1]+g1[n+1]+2*sum)
s=@sprintf("梯形積分計算雙重積分結果 T%d=%.6lf\n",n,ts)
println(s)
s=@sprintf("實際值=%.6lf\n",(1/3))
println(s)

s=@sprintf("誤差值=%.6lf\n",(abs(1/3)- ts))
println(s)

輸出畫面

梯形積分計算雙重積分結果 T10=0.331650

實際值=0.333333

誤差值=0.001683

Julia語言例題4-4使用求 f(x)= (1/x) a=1.0 , b=2.0 求n=?? 辛普森積分法 , 梯形積分法需要幾次

Julia語言例題4-4使用
求 f(x)= (1/x) a=1.0 , b=2.0 求n=??

辛普森積分法需要 n 次 n=8
梯形積分法需要 n 次 n=41

x2
為了方便，符號簡化: ∫ dx = ∫ , h = x2-x0
x0
Thm:

Given [x0,x2] , x1 = (x0+x2)/2, assume f in C^4[x0,x2]

∫f = (h/3) * ( f(x0) + 4f(x1) + f(x2) ) - (h^5/90)*f^(4)(ξ)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~
這項是標準的辛普森法這項是誤差

#========================================================
/* ex4-4.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F4(x::Float64) #// 四次微分後函數
return (24/x^5)
end

function F1(x::Float64) #// 欲微分函數
return (1/x)
end

function F2(x::Float64) #// 二次微分後函數
return (2/x^3)
end

a=1.0
b=2.0
TOL=0.0001
xa=F4(a)
xb=F4(b)
if xa>xb
k=xa
else
k=xb
end

n1= k* (b-a)^5
n2= n1 / (180*TOL)
n3= sqrt(n2)
n4=sqrt(n3)
n=Int32(round(n4+0.5))
if mod(n,2)!=0
n=n+1
end
println("---------------------------------")
println("辛普森積分法需要",n,"次")

#驗證
println("---------------------------------")
sn=0.0
a=1.0
b=2.0
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)

#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
else
sum1=sum1+F1(x)
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)

#s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(sn-(log(2)-log(1))))
println(s)
end
end
print("辛普森積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,sn ,abs(sn-(log(2)-log(1))))
println(s)

#梯形法

a=1.0
b=2.0
TOL=0.0001
xa=F2(a)
xb=F2(b)
if xa>xb
k=xa
else
k=xb
end

n1= k* (b-a)^5
n2= n1 / (12*TOL)
n4= sqrt(n2)
n=Int32(round(n4+0.5))
println("---------------------------------")
println("梯形積分法需要",n,"次")

println("---------------------------------")
sn=0.0
a=1.0
b=2.0
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
result=result+F1(abs(x))
tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f err=%0.6f ",i,x,F1(abs(x)), abs(tn-(log(2)-log(1))))
println(s)

end

#tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("S%d=%lf 誤差=%0.6f\n",n,tn ,abs(tn-(log(2)-log(1))))
println(s)

輸出畫面

---------------------------------
辛普森積分法需要8次
---------------------------------
i= 1 , x=1.13 ---- f(x)=0.888889 err=0.482499 
i= 2 , x=1.25 ---- f(x)=0.800000 err=0.415832 
i= 3 , x=1.38 ---- f(x)=0.727273 err=0.294620 
i= 4 , x=1.50 ---- f(x)=0.666667 err=0.239065 
i= 5 , x=1.63 ---- f(x)=0.615385 err=0.136501 
i= 6 , x=1.75 ---- f(x)=0.571429 err=0.088882 
i= 7 , x=1.88 ---- f(x)=0.533333 err=0.000007 
辛普森積分法  S8=0.693155 誤差=0.000007

---------------------------------
梯形積分法需要41次
---------------------------------
i= 1 , x=1.02 ---- f(x)=0.976190 err=0.651045 
i= 2 , x=1.05 ---- f(x)=0.953488 err=0.627789 
i= 3 , x=1.07 ---- f(x)=0.931818 err=0.605062 
i= 4 , x=1.10 ---- f(x)=0.911111 err=0.582840 
i= 5 , x=1.12 ---- f(x)=0.891304 err=0.561101 
i= 6 , x=1.15 ---- f(x)=0.872340 err=0.539824 
i= 7 , x=1.17 ---- f(x)=0.854167 err=0.518991 
i= 8 , x=1.20 ---- f(x)=0.836735 err=0.498582 
i= 9 , x=1.22 ---- f(x)=0.820000 err=0.478582 
i=10 , x=1.24 ---- f(x)=0.803922 err=0.458975 
i=11 , x=1.27 ---- f(x)=0.788462 err=0.439744 
i=12 , x=1.29 ---- f(x)=0.773585 err=0.420876 
i=13 , x=1.32 ---- f(x)=0.759259 err=0.402357 
i=14 , x=1.34 ---- f(x)=0.745455 err=0.384176 
i=15 , x=1.37 ---- f(x)=0.732143 err=0.366318 
i=16 , x=1.39 ---- f(x)=0.719298 err=0.348775 
i=17 , x=1.41 ---- f(x)=0.706897 err=0.331533 
i=18 , x=1.44 ---- f(x)=0.694915 err=0.314584 
i=19 , x=1.46 ---- f(x)=0.683333 err=0.297917 
i=20 , x=1.49 ---- f(x)=0.672131 err=0.281524 
i=21 , x=1.51 ---- f(x)=0.661290 err=0.265395 
i=22 , x=1.54 ---- f(x)=0.650794 err=0.249522 
i=23 , x=1.56 ---- f(x)=0.640625 err=0.233897 
i=24 , x=1.59 ---- f(x)=0.630769 err=0.218512 
i=25 , x=1.61 ---- f(x)=0.621212 err=0.203361 
i=26 , x=1.63 ---- f(x)=0.611940 err=0.188435 
i=27 , x=1.66 ---- f(x)=0.602941 err=0.173729 
i=28 , x=1.68 ---- f(x)=0.594203 err=0.159237 
i=29 , x=1.71 ---- f(x)=0.585714 err=0.144951 
i=30 , x=1.73 ---- f(x)=0.577465 err=0.130867 
i=31 , x=1.76 ---- f(x)=0.569444 err=0.116978 
i=32 , x=1.78 ---- f(x)=0.561644 err=0.103279 
i=33 , x=1.80 ---- f(x)=0.554054 err=0.089765 
i=34 , x=1.83 ---- f(x)=0.546667 err=0.076432 
i=35 , x=1.85 ---- f(x)=0.539474 err=0.063274 
i=36 , x=1.88 ---- f(x)=0.532468 err=0.050287 
i=37 , x=1.90 ---- f(x)=0.525641 err=0.037467 
i=38 , x=1.93 ---- f(x)=0.518987 err=0.024809 
i=39 , x=1.95 ---- f(x)=0.512500 err=0.012309 
i=40 , x=1.98 ---- f(x)=0.506173 err=0.000037 
梯形積分法  S41=0.693184 誤差=0.000037

Julia語言例題4-4利用辛普森理則（Simpson's Rule)與梯形法計算定積分

Julia語言例題4-4利用辛普森理則（Simpson's Rule) 與梯形法
# (a) 計算 exp( 1/x) 在[1, 2]的定積分
# (b) 計算 ( 1/x) 在[1, 2]的定積分

#========================================================
/* ex4-4.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F1(x::Float64) #// 欲微分函數
return (exp(1/x))
end

function F2(x::Float64) #// 欲微分函數
return (1/x)
end

a=1.0
b=2.0
n=4
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F1(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
println(s)
else
sum1=sum1+F1(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F1(x))
println(s)
end
end
sn= (h/3.0)*(F1(a)+F1(b)+2.0*sum2+4.0*sum1)
print("辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

println("---------------------------------")
n=10
a=1.0
b=2.0
n=4
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f",i,x,F1(abs(x)))
println(s)
result=result+F1(abs(x))
end

tn=(h/2.0)*(F1(abs(a))+F1(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

println("---------------------------------")
n=10
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F2(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F2(x))
println(s)
else
sum1=sum1+F2(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F2(x))
println(s)
end
end
sn= (h/3.0)*(F2(a)+F2(b)+2.0*sum2+4.0*sum1)
print("辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

println("---------------------------------")
n=10
a=1.0
b=2.0
n=10
h=(b-a)/n
x=a
result=0.0
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.6f",i,x,F2(abs(x)))
println(s)
result=result+F2(abs(x))
end

tn=(h/2.0)*(F2(abs(a))+F2(abs(b))+2.0*result)

print("梯形積分法 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

輸出畫面

i= 1 ,x=1.250 --- F(x)=2.225541
i= 2 ,x=1.500 --- F(x)=1.947734
i= 3 ,x=1.750 --- F(x)=1.770795
辛普森積分法  S4=2.020651

---------------------------------
i= 1 , x=1.25 ---- f(x)=2.225541
i= 2 , x=1.50 ---- f(x)=1.947734
i= 3 , x=1.75 ---- f(x)=1.770795
梯形積分法  T4=2.0318929

---------------------------------
i= 1 ,x=1.100 --- F(x)=0.909091
i= 2 ,x=1.200 --- F(x)=0.833333
i= 3 ,x=1.300 --- F(x)=0.769231
i= 4 ,x=1.400 --- F(x)=0.714286
i= 5 ,x=1.500 --- F(x)=0.666667
i= 6 ,x=1.600 --- F(x)=0.625000
i= 7 ,x=1.700 --- F(x)=0.588235
i= 8 ,x=1.800 --- F(x)=0.555556
i= 9 ,x=1.900 --- F(x)=0.526316
辛普森積分法  S10=0.693150

---------------------------------
i= 1 , x=1.10 ---- f(x)=0.909091
i= 2 , x=1.20 ---- f(x)=0.833333
i= 3 , x=1.30 ---- f(x)=0.769231
i= 4 , x=1.40 ---- f(x)=0.714286
i= 5 , x=1.50 ---- f(x)=0.666667
i= 6 , x=1.60 ---- f(x)=0.625000
i= 7 , x=1.70 ---- f(x)=0.588235
i= 8 , x=1.80 ---- f(x)=0.555556
i= 9 , x=1.90 ---- f(x)=0.526316
梯形積分法  T10=0.6937714

Julia語言例題4-5 利用辛普森理則（Simpson's Rule) 計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

Julia語言例題4-5 利用辛普森理則（Simpson's Rule)
#計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

#========================================================
/* ex4-5.jl based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
========================================================#
using Printf

function F(x::Float64) #// 欲微分函數
return (( (1+x^3)^ (1/3) )^2)

end
a=0.0
b=2.0
n=10
m=n/2
h=(b-a)/n
sum1=0.0
sum2=0.0
for i=1:2*m-1
x=a+i*h;
if(i%2==0)
sum2=sum2+F(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F(x))
println(s)
else
sum1=sum1+F(x)
s=@sprintf("i=%2d ,x=%0.3f --- F(x)=%0.6f",i,x,F(x))
println(s)
end
end
sn= pi *(h/3.0)*(F(a)+F(b)+2.0*sum2+4.0*sum1)
print("\n\n辛普森積分法 ")
s=@sprintf("S%d=%lf\n",n,sn)
println(s)

輸出畫面

i= 1 ,x=0.200 --- F(x)=1.005326
i= 2 ,x=0.400 --- F(x)=1.042224
i= 3 ,x=0.600 --- F(x)=1.139259
i= 4 ,x=0.800 --- F(x)=1.317350
i= 5 ,x=1.000 --- F(x)=1.587401
i= 6 ,x=1.200 --- F(x)=1.952374
i= 7 ,x=1.400 --- F(x)=2.411148
i= 8 ,x=1.600 --- F(x)=2.961326
i= 9 ,x=1.800 --- F(x)=3.600520


辛普森積分法  S10=12.325078

Julia語言例題EX4-2 請用梯形法求統計學常態分布的面積

f(x) = \frac{1}{\sqrt{2\pi}} \, \exp\left(-\frac{x^2}{2} \right)

。

a=-3.5 b=3.0 n=10 的積分面積

程式
#================================================
* ex4-2.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

=================================================#

using Printf

function f(x::Float64) #// 欲微分函數
return return (1.0/exp(x*x/2))
end

result=0.0
n=10
a=-3.5
b=-3.0

h=(b-a)/n
x=a
for i=1:n-1
x=x+h
s=@sprintf("i=%2d , x=%0.2f ---- f(x)=%0.4f",i,x,f(abs(x)))
println(s)
result=result+f(abs(x))
end

tn=(1.0/sqrt(2* pi ))*(h/2.0)*(f(abs(a))+f(abs(b))+2.0*result)

println("\nf(x)=(1.0/exp(x*x/2)) a=0.0 ,b=1.0 n=10 ")
s=@sprintf("T%d=%0.7lf\n",n,tn)
println(s)

ET1= (b^2 -1)* exp(- (b^2)/2)
#println(ET1)
ET= (ET1 * (b-a)^3) / (12*n^2)
s=@sprintf("梯形法的誤差範圍= %0.7f",ET)
println(s)

輸出畫面

i= 1 , x=-3.45 ---- f(x)=0.0026
i= 2 , x=-3.40 ---- f(x)=0.0031
i= 3 , x=-3.35 ---- f(x)=0.0037
i= 4 , x=-3.30 ---- f(x)=0.0043
i= 5 , x=-3.25 ---- f(x)=0.0051
i= 6 , x=-3.20 ---- f(x)=0.0060
i= 7 , x=-3.15 ---- f(x)=0.0070
i= 8 , x=-3.10 ---- f(x)=0.0082
i= 9 , x=-3.05 ---- f(x)=0.0095

f(x)=(1.0/exp(x*x/2)) a=0.0 ,b=1.0 n=10 
T10=0.0011194

梯形法的誤差範圍= 0.0000093

Julia語言利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

Julia語言例題4-1利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

#================================================
* ex4-1.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

=================================================#

using Printf

function f(x::Float64) #// 欲微分函數
return exp(x)
end

result=0.0
n=10
a=0.0
b=1.0

h=(b-a)/n
x=a
for i=1:n-1
x=x+h;
s=@sprintf("x=%0.1f ---- f(x)=%0.4f",x,f(abs(x)))
println(s)
result=result+f(abs(x))
end

tn=(h/2.0)*( f(abs(a)) + f(abs(b)) + 2.0*result )

println("\nf(x)=exp(x) a=0.0 ,b=1.0 n=10 ")
s=@sprintf("T%d=%10.6lf\n",n,tn)
println(s)

s=@sprintf("梯形積分與實際值的誤差值=%0.5f",abs(tn-(f(1.0)-f(0.0))))
println(s)

ET= (f(b)* (b-a)^3) / (12*n^2)
s=@sprintf("梯形法的誤差範圍= %0.5f",ET)
println(s)

輸出畫面

x=0.1 ---- f(x)=1.1052
x=0.2 ---- f(x)=1.2214
x=0.3 ---- f(x)=1.3499
x=0.4 ---- f(x)=1.4918
x=0.5 ---- f(x)=1.6487
x=0.6 ---- f(x)=1.8221
x=0.7 ---- f(x)=2.0138
x=0.8 ---- f(x)=2.2255
x=0.9 ---- f(x)=2.4596

f(x)=exp(x) a=0.0 ,b=1.0 n=10 
T10=  1.719713

梯形積分與實際值的誤差值=0.00143

梯形法的誤差範圍= 0.00227

Julia語言例題4-3 梯形法求標準常態分配的面積

f(x) = {1 \over \sigma\sqrt{2\pi} }\,e^{- {{(x-\mu )^2 \over 2\sigma^2}}}

f(x) = \frac{1}{\sqrt{2\pi}} \, \exp\left(-\frac{x^2}{2} \right)

a=-5 , b=5 , n=100

#================================================
* ex4-3.jl based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*
梯形法 T'= (h/2) [ f(a) + f(b) + 2 sigma i=1 to n-1 f(xi)]

=================================================#

using Printf

function f(x::Float64) #// 欲微分函數
return (1.0/exp(x*x/2))
end

result=0.0
n=100
a=-5.0
b=5.0

h=(b-a)/n
x=a
for i=1:n
x=x+h;
result=result+f(abs(x))
end

tn=(1.0/sqrt(2* pi ))*(h/2.0)*(f(abs(a))+f(abs(b))+2.0*result)

s=@sprintf("T%d=%10.6lf\n",n,tn)
println(s)

輸出畫面

T100=  1.000000

Julia語言習題3-4請使用中央差近似法計算一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

Julia語言習題3-4 觀察下列數據
x = 1.00 , 1.01 , 1.02
f(x)=1.27 , 1.32 , 1.38
請使用中央差近似法計算
一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

程式
#中央近似差方法
# f' i = (f(i+1) - f(i-1)) / 2h
# f'' i = ( fi+1 - 2 fi + f i-1) / h^2
#============================================
x f(x) f'(x) f''(x)
1.00 1.27
1.005
1.01 1.32
1.015
1.02 1.38
============================================#

using Printf

x=[ 1.00 , 1.01 , 1.02 ]
f= [ [1.27 ,0.0 , 0.0 ],
[1.32 ,0.0 , 0.0 ],
[1.38 ,0.0 , 0.0 ] ]

xa=[1.005 , 1.015 ,1.01]
h=abs(x[2]-x[1])
i=1
s=@sprintf("f(%0.3f)的微分.......",xa[i])
println(s)
fd1=(f[i+1][1]- f[i][1]) / (h)
s=@sprintf("一次微分值= %0.4f ",fd1)
f[i][2]=fd1
println(s)

println("\n")
i=2
s=@sprintf("f(%0.3f)的微分.......",xa[i])
println(s)

fd1=(f[i+1][1]- f[i][1]) / (h)
s=@sprintf("一次微分值= %0.4f ",fd1)
f[i][2]=fd1
println(s)

println("\n")
i=2
s=@sprintf("f(%0.3f)的二次微分.......",xa[i+1])
println(s)

fd2=(f[i+1][1]- 2*f[i][1] + f[i-1][1]) / (h*h)
s=@sprintf("二次微分值= %0.4f ",fd2)
f[i][3]=fd2
println(s)

輸出畫面

f(1.005)的微分.......
一次微分值= 5.0000 


f(1.015)的微分.......
一次微分值= 6.0000 


f(1.010)的二次微分.......
二次微分值= 100.0000

Julia語言習題3-3 等距的函數f(x)的微分近似值

Julia語言習題3-3 等距的函數f(x)的微分近似值

先使用牛頓向前的內插多項式
x= f(x)= f ' (x)= f ''(x)=
=============================================
0.15 0.1761 2.8775 -17 前向
0.17 0.2304 2.5675 -14.75 中央公式相同
0.19 0.2788 2.295 -12.5 中央公式相同
0.21 0.3222 2.06 -10.25 後向
h=0.02

程式

using Printf

x=[ 0.15 , 0.17 , 0.19 , 0.21 ]

f= [ [0.1761 ,0.0 , 0.0 ],

[0.2304 ,0.0 , 0.0 ],

[0.2788 ,0.0 , 0.0 ],

[0.3222 ,0.0 , 0.0 ] ]

h=abs(x[1]-x[2])

i=1

s=@sprintf("f(%0.3f)的微分.......",x[i])

println(s)

fd1=( 2*f[i+3][1] - 9*f[i+2][1] + 18*f[i+1][1] - 11*f[i][1] ) / (6*h)

fd2=( - f[i+3][1] + 4*f[i+2][1] + -5*f[i+1][1] + 2*f[i][1] ) / (h*h)

s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)

f[i][2]=fd1

f[i][3]=fd2

println(s)

println("\n")

i=2

s=@sprintf("f(%0.3f)的微分.......",x[i])

println(s)

fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)

fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)

s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)

f[i][2]=fd1

f[i][3]=fd2

println(s)

println("\n")

i=3

s=@sprintf("f(%0.3f)的微分.......",x[i])

println(s)

fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)

fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)

s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)

f[i][2]=fd1

f[i][3]=fd2

println(s)

println("\n")

i=4

s=@sprintf("f(%0.3f)的微分.......",x[i])

println(s)

fd1=( 11*f[i][1] - 18*f[i-1][1] + 9*f[i-2][1] - 2*f[i-3][1] ) / (6*h)

fd2=( 2*f[i][1] - 5*f[i-1][1] + 4*f[i-2][1] - f[i-3][1] ) / (h*h)

s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)

f[i][2]=fd1

f[i][3]=fd2

println(s)

println("\n")

println("x\tf(x)\t\tf'(x)\t\tf''(x)")

for i=1:4

print(x[i],"\t")

s=@sprintf("%0.5f\t\t%0.5f\t\t%0.5f",f[i][1],f[i][2],f[i][3])

println(s)

end

輸出結果

f(0.150)的微分.......
一次微分值= 2.8775 , 二次微分值= -17.0000 


f(0.170)的微分.......
一次微分值= 2.5675 , 二次微分值= -14.7500 


f(0.190)的微分.......
一次微分值= 2.2950 , 二次微分值= -12.5000 


f(0.210)的微分.......
一次微分值= 2.0600 , 二次微分值= -10.2500 


x f(x)  f'(x)  f''(x)
0.15 0.17610  2.87750  -17.00000
0.17 0.23040  2.56750  -14.75000
0.19 0.27880  2.29500  -12.50000
0.21 0.32220  2.06000  -10.25000

2019年2月28日星期四

Julia語言習題3-6請選擇適當的微分近似法完成下表：

Julia語言習題3-6請選擇適當的微分近似法完成下表：
x = 0 ,2 ,3 ,5
f(x)= 1 ,7 , 25 ,121
f'(x) = ? ? ? ?

程式
using Printf

x=[ 0 , 2 , 3 , 5]
f= [ [1.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[7.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[25.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[121.0 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ] ]

n=length(x)
println(" Divided Difference Table: ")
println("=============================")
for j=2:n
for i=1:n-j+1
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
end
end

print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2), ......................\n")
for i=1:n
s=@sprintf("%d\t%8.5f",i,x[i])
print(s)
for j=1:n-i+1
s=@sprintf("\t%8.5f",f[i][j])
print(s)
end
println()

end

d=[0.0 for i=1:n-1]
for i=1:n-1
d[i]=f[1][i+1]
end

s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)

xa=[ 0.0 , 2.0 , 3.0 , 5.0]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
pn[i]=d[1]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
for j=1:2
#println(xb," ", x[j])
dx=dx+(xb-x[j])
end
#println(dx)
dx=d[2]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
dx=d[3]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

輸出畫面

Divided Difference Table: 
=============================
i x(i)  f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  0.00000  1.00000  3.00000  5.00000  1.00000
2  2.00000  7.00000 18.00000 10.00000
3  3.00000 25.00000 48.00000
4  5.00000 121.00000

牛頓一階微分前向差除表[3.0, 5.0, 1.0]

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = 3.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 3.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 3.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 3.00000

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = -7.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 13.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 23.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 43.00000

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.000 , P'n( 0.000 ) = -1.00000
========================================================
xa[2]= 2.000 , P'n( 2.000 ) = 11.00000
========================================================
xa[3]= 3.000 , P'n( 3.000 ) = 26.00000
========================================================
xa[4]= 5.000 , P'n( 5.000 ) = 74.00000

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)
微分的近似值完成表中的 f'(x) 與 f''(x)
x= f(x)= f'(x)= f''(x)=
=============================
0.15 0.1761 ?? ??
0.17 0.2304 ?? ??
0.19 0.2788 ?? ??
0.21 0.3222 ?? ??
=============================
h=0.02

using Printf

x=[ 0.15 , 0.17 , 0.19 , 0.21 ]
f= [ [0.1761 ,0.0 , 0.0 ],
[0.2304 ,0.0 , 0.0 ],
[0.2788 ,0.0 , 0.0 ],
[0.3222 ,0.0 , 0.0 ] ]

h=abs(x[1]-x[2])
i=1
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=( 2*f[i+3][1] - 9*f[i+2][1] + 18*f[i+1][1] - 11*f[i][1] ) / (6*h)
fd2=( - f[i+3][1] + 4*f[i+2][1] + -5*f[i+1][1] + 2*f[i][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=2
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)
fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=3
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)
fd1=(f[i+1][1]- f[i+-1][1]) / (2*h)
fd2=(f[i+1][1] - 2*f[i][1] + f[i-1][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
i=4
s=@sprintf("f(%0.3f)的微分.......",x[i])
println(s)

fd1=( 11*f[i][1] - 18*f[i-1][1] + 9*f[i-2][1] - 2*f[i-3][1] ) / (6*h)
fd2=( 2*f[i][1] - 5*f[i-1][1] + 4*f[i-2][1] - f[i-3][1] ) / (h*h)
s=@sprintf("一次微分值= %0.4f , 二次微分值= %0.4f ",fd1,fd2)
f[i][2]=fd1
f[i][3]=fd2
println(s)

println("\n")
println("x\tf(x)\t\tf'(x)\t\tf''(x)")
for i=1:4
print(x[i],"\t")
s=@sprintf("%0.5f\t\t%0.5f\t\t%0.5f",f[i][1],f[i][2],f[i][3])
println(s)
end

輸出畫面

f(0.150)的微分.......
一次微分值= 2.8775 , 二次微分值= -17.0000 


f(0.170)的微分.......
一次微分值= 2.5675 , 二次微分值= -14.7500 


f(0.190)的微分.......
一次微分值= 2.2950 , 二次微分值= -12.5000 


f(0.210)的微分.......
一次微分值= 2.0600 , 二次微分值= -10.2500

x     f(x)     f'(x)     f''(x)

0.15 0.17610  2.87750  -17.00000
0.17 0.23040  2.56750  -14.75000
0.19 0.27880  2.29500  -12.50000
0.21 0.32220  2.06000  -10.25000

Julia語言例題3-3 不等距的函數f(x)的微分近似值先使用牛頓向前的內插多項式

Julia語言例題3-3 不等距的函數f(x)的微分近似值
先使用牛頓向前的內插多項式

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
x f(x) 第一項P'n(x) 前二項P'n(x)　　前三項P'n(x)
0.5 0.4794
0.6 0.5646
0.8 0.7174
1.05 0.8674
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝====
P'(n) 第1項 f0,1

P'(n) 第2項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )

P'(n) 第3項 f0,1+f0,1,2 *( (x-x[1]) + (x-x[2]) )
+ f0,1,2,3 * ( (x-x[1])*(x-x[2]) ) + ( (x-x[2])*(x-x[3]) ) + ( (x-x[3])*(x-x[1]) )

using Printf

x=[ 0.5 , 0.6 , 0.8 , 1.05]
f= [ [0.4794 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.5646 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.7174 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ],
[0.8674 ,0.0 , 0.0 ,0.0 ,0.0 ,0.0 ] ]

n=length(x)
println(" Divided Difference Table: ")
println("=============================")
for j=2:n
for i=1:n-j+1
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j-1]-x[i])
end
end

print("i\tx(i)\t\tf(i)\t\tf(i,i+1)\tf(i,i+1.i+2), ......................\n")
for i=1:n
s=@sprintf("%d\t%8.5f",i,x[i])
print(s)
for j=1:n-i+1
s=@sprintf("\t%8.5f",f[i][j])
print(s)
end
println()

end

d=[0.0 for i=1:n-1]
for i=1:n-1
d[i]=f[1][i+1]
end

s=@sprintf("\n牛頓一階微分前向差除表")
println(s,d)

xa=[ 0.5 , 0.6 , 0.8 , 1.05]
pn=[0.0 , 0.0 , 0.0 ,0.0]
s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
pn[i]=d[1]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
for j=1:2
#println(xb," ", x[j])
dx=dx+(xb-x[j])
end
#println(dx)
dx=d[2]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

s=@sprintf("\nTHE RESULTS OF INTERPOLATION:\n")
print(s)
for i=1:length(xa)
xb=xa[i]
dx=0.0
dx=(xb-x[1])*(xb-x[2])+ (xb-x[2])*(xb-x[3])+ (xb-x[3])*(xb-x[1])
dx=d[3]*dx
pn[i]=pn[i]+dx
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , P'n( %0.3f ) = %0.5f",i,xa[i],xa[i],pn[i] )
println(s)
end

println("\n\n")
for i=1:length(xa)
xb=xa[i]
println("========================================================")
s=@sprintf("xa[%1d]= %0.3f , 誤差 f'(x)-P'(n) = %0.5f",i,xa[i],abs(cos(xb)-pn[i]) )
println(s)
end

輸出畫面

i f(i)  f(i,i+1) f(i,i+1.i+2),  ......................
1  0.50000  0.47940  0.85200 -0.29333 -0.12929
2  0.60000  0.56460  0.76400 -0.36444
3  0.80000  0.71740  0.60000
4  1.05000  0.86740

牛頓一階微分前向差除表[0.852, -0.293333, -0.129293]

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.85200
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.85200
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.85200
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.85200

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.88133
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82267
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.70533
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.55867

THE RESULTS OF INTERPOLATION:
========================================================
xa[1]= 0.500 , P'n( 0.500 ) = 0.87745
========================================================
xa[2]= 0.600 , P'n( 0.600 ) = 0.82525
========================================================
xa[3]= 0.800 , P'n( 0.800 ) = 0.69758
========================================================
xa[4]= 1.050 , P'n( 1.050 ) = 0.49434

========================================================
xa[1]= 0.500 , 誤差 f'(x)-P'(n) = 0.00013
========================================================
xa[2]= 0.600 , 誤差 f'(x)-P'(n) = 0.00008
========================================================
xa[3]= 0.800 , 誤差 f'(x)-P'(n) = 0.00087
========================================================
xa[4]= 1.050 , 誤差 f'(x)-P'(n) = 0.00323

2019年2月25日星期一

Julia語言利用3點微分方式求f'(x)

Julia語言利用3點微分方式求f'(x)
3點向前
fd1=((-f[i+2]+ 4*f[i+1]-3*f[i])/ (2*h))
3點向後
fd2=((3*f[m] -4*f[m-1]+ f[m-2]) / (2*h))

x f(x) f'(x)
======================
-0.3 -0.20431 ???
-0.1 -0.08993 ???
0.1 0.11007 ???
0.3 0.39569 ???
======================

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:

sin cos tan cot sec csc
sinh cosh tanh coth sech csch
asin acos atan acot asec acsc
asinh acosh atanh acoth asech acsch
sinc cosc atan2
=====================================================#

x=[-0.3 , -0.1, 0.1 , 0.3 ]
f=[-0.20431 , -0.08993 , 0.11007 , 0.39569]

println("3點向前差")
for i=1: (length(x)-2)
h= abs(x[i+1]-x[i])
fd1=((-f[i+2]+ 4*f[i+1]-3*f[i])/ (2*h))
s=@sprintf("3點向前差 Three-point fordward Diff : fd(%0.5lf) =%0.5lf \n",x[i], fd1)
println(s)
end

println("3點向後差")
for m=length(x):-1:3
h= abs(x[m]-x[m-1])
fd2=((3*f[m] -4*f[m-1]+ f[m-2]) / (2*h))
s=@sprintf("3點向後差 Three-point bachfoward Diff : fd(%0.5lf) =%0.5lf \n",x[m], fd2)
println(s)
end

輸出畫面

3點向前差
3點向前差 Three-point fordward Diff : fd(-0.30000) =0.35785 

3點向前差 Three-point fordward Diff : fd(-0.10000) =0.78595 

3點向後差
3點向後差 Three-point bachfoward Diff : fd(0.30000) =1.64215 

3點向後差 Three-point bachfoward Diff : fd(0.10000) =1.21405

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值
h=0.001 , 0.005 , 0.01 , 0.05 , 0.1 ,0.5

def func( x): # // 欲微分函數
return math.sin(x)

def FordDiff( x, h, fx):
# // 前差微分
return ( ( fx(x+h) - fx(x) ) / h)

def BackDiff(x, h, fx):
# // 後差微分
return (( fx(x) - fx(x-h) ) / h)

def MidDiff(x, h, fx):
# // 中差微分
return (0.5 * ( fx(x+h) - fx(x-h) ) / h)

using Printf

#=====================================================
Trigonometric and hyperbolic functions
All the standard trigonometric and hyperbolic functions are also defined:

sin cos tan cot sec csc
sinh cosh tanh coth sech csch
asin acos atan acot asec acsc
asinh acosh atanh acoth asech acsch
sinc cosc atan2
=====================================================#

function func(x::Float64) #// 欲微分函數
return sin(x)
end

function FordDiff(x::Float64 , h::Float64, fx::typeof(func)) #//前差微分
a=fx(x+h)
b=fx(x)
return ((a-b) / h)
end

function BackDiff(x::Float64 , h::Float64, fx::typeof(func)) #//後差微分
return ((fx(x) - fx(x-h) ) / h )
end

function MidDiff(x::Float64 , h::Float64, fx::typeof(func)) #//中差微分
return (0.5 * ( fx(x+h) - fx(x-h) ) / h)
end

#=====================================================#
hx=[0.001 , 0.005 , 0.01 , 0.05 , 0.1 , 0.5]

for m=1:length(hx)
x=1.0
h=hx[m]
println("h=",h,"\n")
answer = cos(x) # // 答案
cal = FordDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("FordDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)

cal = BackDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("BackDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)

cal = MidDiff(x, h, func)
delta = (cal - answer)/answer
s=@sprintf("MidDiff : %0.5lf, delta = %0.5lf %%\n",cal,delta)
println(s)
println("================================================")
end

輸出畫面

h=0.001

FordDiff : 0.53988, delta = -0.00078 %

BackDiff : 0.54072, delta = 0.00078 %

MidDiff : 0.54030, delta = -0.00000 %

================================================
h=0.005

FordDiff : 0.53820, delta = -0.00390 %

BackDiff : 0.54240, delta = 0.00389 %

MidDiff : 0.54030, delta = -0.00000 %

================================================
h=0.01

FordDiff : 0.53609, delta = -0.00780 %

BackDiff : 0.54450, delta = 0.00777 %

MidDiff : 0.54029, delta = -0.00002 %

================================================
h=0.05

FordDiff : 0.51904, delta = -0.03934 %

BackDiff : 0.56111, delta = 0.03851 %

MidDiff : 0.54008, delta = -0.00042 %

================================================
h=0.1

FordDiff : 0.49736, delta = -0.07947 %

BackDiff : 0.58144, delta = 0.07614 %

MidDiff : 0.53940, delta = -0.00167 %

================================================
h=0.5

FordDiff : 0.31205, delta = -0.42246 %

BackDiff : 0.72409, delta = 0.34016 %

MidDiff : 0.51807, delta = -0.04115 %

================================================

訂閱：文章 (Atom)

alex9ufo 聰明人求知心切

2019年3月7日星期四

Julia語言例題1-1 若已知下面二點請問 x=1.5時則f(1.5)= ??

2019年3月5日星期二

Julia語言例題1-5 已知函數f(x)=ln(x)請使用Lagrange 內插法求f(1.5) ,f(2.5) f(3.5)的P(x)與f(x) 分別用x=1.5 ,2.5 ,3.5 計算f(x)的e(x)誤差值

2019年3月3日星期日

Julia語言例題EX4-7-B利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 n-4 並比較誤差值

2019年3月2日星期六

Julia語言例題EX4-7-A利用梯形法求雙重積分 f (x,y) = x* exp(y) a=0.0 , b=1.0 , c(x) =0.0 d(x)= x 取 n=2 與 n=4 二者並比較誤差值

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=10 求積分值

Julia語言例題4-6 雙重定積分 a=0.0 , b=1.0 , c=1.0 ,d=2.0 , f(x,y)= x^2 y 取n=2 求積分值

2019年3月1日星期五

Julia語言範例EX4-8 梯行積分法雙重積分計算

Julia語言例題4-4使用求 f(x)= (1/x) a=1.0 , b=2.0 求n=?? 辛普森積分法 , 梯形積分法需要幾次

Julia語言例題4-4利用辛普森理則（Simpson's Rule)與梯形法計算定積分

Julia語言例題4-5 利用辛普森理則（Simpson's Rule) 計算 x=0 to 2 , 以曲線 y= (1+x^3) ^(1/3) 沿 z 軸旋轉一周的體積

Julia語言例題EX4-2 請用梯形法求統計學常態分布的面積

Julia語言利用梯形積分法求f (x) = exp(x) a=0.0 , b=1.0 n=10的積分值

Julia語言例題4-3 梯形法求標準常態分配的面積

Julia語言習題3-4請使用中央差近似法計算一次微分 f '(1.005) , f ' (1.015) 與二次微分 f '' (1.01) 之值

Julia語言習題3-3 等距的函數f(x)的微分近似值

2019年2月28日星期四

Julia語言習題3-6請選擇適當的微分近似法完成下表：

Julia語言微分的近似值習題3-3 完成表中的 f'(x) 與 f''(x)

Julia語言例題3-3 不等距的函數f(x)的微分近似值先使用牛頓向前的內插多項式

2019年2月25日星期一

Julia語言利用3點微分方式求f'(x)

Julia語言例題3-1 利用向前差　向後差　中央差　求f(x)=sin(x)的值

MQTT 協定與 Modbus 通訊的遠端監控與控制系統架構

2019年3月7日 星期四

2019年3月5日 星期二

2019年3月3日 星期日

2019年3月2日 星期六

2019年3月1日 星期五

2019年2月28日 星期四

2019年2月25日 星期一

2019年3月7日星期四

2019年3月5日星期二

2019年3月3日星期日

2019年3月2日星期六

2019年3月1日星期五

2019年2月28日星期四

2019年2月25日星期一