微分方程数值解作业 2
微分方程
Problem 1
To solve (2.4)
We can always obtain the equation at non-boundary positions,
The above FDE approximation has an error of
Question (a)
The first condition
and in Eqn. becomes
For the second condition
substitute with an one-sided approximation
we have
And the final equation becomes
The whole FDE system in matrix form
Question (b)
Therefore
FDE system in matrix form
Problem 2
Question (a)
Use central difference to approximate with
Substitute back
Question (b)
Given the differerntial equation
Simply susbtitute with in the above equation, and discard the higher order terms, we have
Question (c)
When
we will be able to seperate from , and the difference equation becomes
Collect terms on the left hand side
The above equation is linear, and could be written in matrix form. Denote
Assume the boundary conditions are and , we can obtain Equation at and
Finally, the FDE system is expressed in explicit matrix form
Question (d)
Therefore, when satisfies , Equation (2.4) will not contain term, and can be solve by the method in Question (b). That requires to be a constant.
Problem 3
At , use a third-ordered forward difference to approximiate and
And use a forth-ordered forward difference to approximiate
Substitute back to the equation we have one equation
with 5 variables to .
At to , use a forth-ordered central difference to approximiate and
and a forth-ordered central difference to approximiate
substitute them back we have equations
with variables from to .
Also, the boundry conditions provide 3 equations
Put Equations together, we could obtain a second-order non-linear FDE system with equations and variables to .
Problem 4
Question (a)
1 2 3 4 eqn = \[Epsilon] D [ y [ x ] , { x , 2 } ] - D [ y [ x ] , x ] == - 1 && y [ 0 ] == 1 && y [ 1 ] == 3 ; sol = ( 1 + # + ( E ^ ( # / \[Epsilon] ) - 1 ) / ( E ^ ( 1 / \[Epsilon] ) - 1 ) ) &; eqn /. { y -> sol }
Question (b)
Collect terms
At boundries
FDE system in matrix form
Since we have used first-ordered backward difference to approximate , the truncation error is .
Question (c)
It is obvious that the matrix in Equation meets the second condition of Theorem 2.1, therefore the matrix is invertible.
If we use Equation (2.16) to build the FDE system, Theorem 2.2 guarantees the existence of a unique solution if . That is, .
Question (d)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 tridiagonalSolve [ b_ , a_ , c_ , z_ ] := Module [ { y , v , w , n } , n = Length [ a ] ; w = a [ [ 1 ] ] ; y = ConstantArray [ 0 , n ] ; y [ [ 1 ] ] = z [ [ 1 ] ] / w ; v = ConstantArray [ 0 , n ] ; Do [ v [ [ i ] ] = c [ [ i - 1 ] ] / w ; w = a [ [ i ] ] - b [ [ i ] ] v [ [ i ] ] ; y [ [ i ] ] = ( z [ [ i ] ] - b [ [ i ] ] y [ [ i - 1 ] ] ) / w ; , { i , 2 , n } ] ; Do [ y [ [ i ] ] = y [ [ i ] ] - v [ [ i + 1 ] ] y [ [ i + 1 ] ] ; , { i , n - 1 , 1 , - 1 } ] ; Return [ y ] ; ] ;
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 solve1 [ \[Epsilon] _ , n_ ] := Module [ { c , a , b , h , x , y , z } , h = 1 / ( n + 1 ) ; x = Range [ h , 1 - h , h ] ; b = ConstantArray [ h + \[Epsilon] , n ] ; a = ConstantArray [ - ( h + 2 \[Epsilon] ) , n ] ; c = ConstantArray [ \[Epsilon] , n ] ; z = ConstantArray [ - h ^ 2 , n ] ; z [ [ 1 ] ] = - h ^ 2 - h - \[Epsilon] ; z [ [ n ] ] = - h ^ 2 - 3 \[Epsilon] ; y = tridiagonalSolve [ b , a , c , z ] ; Return [ Interpolation [ Join [ { { 0 , 1 } , { 1 , 3 } } , Transpose [ { x , y } ] ] , InterpolationOrder -> 1 ] ] ; ] ; solve2 [ \[Epsilon] _ , n_ ] := Module [ { c , a , b , h , x , y , z } , h = 1 / ( n + 1 ) ; x = Range [ h , 1 - h , h ] ; b = ConstantArray [ \[Epsilon] + h / 2 , n ] ; a = ConstantArray [ - 2 \[Epsilon] , n ] ; c = ConstantArray [ \[Epsilon] - h / 2 , n ] ; z = ConstantArray [ - h ^ 2 , n ] ; z [ [ 1 ] ] = - h ^ 2 - b [ [ 1 ] ] ; z [ [ n ] ] = - h ^ 2 - 3 c [ [ n ] ] ; y = tridiagonalSolve [ b , a , c , z ] ; Return [ Interpolation [ Join [ { { 0 , 1 } , { 1 , 3 } } , Transpose [ { x , y } ] ] , InterpolationOrder -> 1 ] ] ; ] ;
1 2 3 4 5 6 7 8 9 10 11 solveAndPlot [ eps_ , n_ ] := Module [ { s1 , s2 } , s1 = solve1 [ eps , n ] ; s2 = solve2 [ eps , n ] ; Plot [ { symSol /. \[Epsilon] -> eps , s1 [ x ] , s2 [ x ] } , { x , 0 , 1 } , PlotLegends -> Placed [ { "Symbolic" , "Method from (b)" , "Equation (2.16)" } , Right ] , PlotLabel -> StringForm [ "Solution for \[Epsilon]=``, N=``" , eps , n ] ] ] ; GraphicsColumn [ { solveAndPlot [ 0.1 , 10 ] , solveAndPlot [ 0.1 , 20 ] , solveAndPlot [ 0.1 , 40 ] } ]
Question (e)
1 2 GraphicsColumn [ { solveAndPlot [ 0.01 , 10 ] , solveAndPlot [ 0.01 , 20 ] , solveAndPlot [ 0.01 , 40 ] } ]
Question (f)
It can be observed that both methods are converging to the symbolic solution as increases. The method in (b) is less accurate than Equation (2.16) when is . However, when gets smaller, (b) is still stable but Equation (2.16) oscillates. When is even smaller, (2.16) oscillates more severely.
The oscillation in (2.16) happens when is violated. Therefore, when is small, the method in (b) better than (2.16) since it converges regardless of and .