Download Numerical Solution of Differential Equations with Finite Difference and Iterative and more Study notes Calculus in PDF only on Docsity!
Math 128A Lab 3: Nonlinear BVPs
Adrian Down
February 27, 2006
Part 1
Discussion
Derivation
I derive a general iterative method to solve BVPs. I consider the case of a differential equation, which I write as
Aφ + g(φ) = 0
A represents the coefficients of the linear components of the differential equa- tion, and g(φ) represents the nonlinear dependence in the differential equa- tion. I assume that I have some proposed solution to the differential equation, call it φk. However, the solution is not exact,
Aφk + g(φk) 6 = 0
I would like to find a small correction, sk, such that
A(φk + sk) + g(φk + sk) = 0
Since A is linear, it distributes over φk and sk. The g term is nonlinear, however. I use a linear approximation for g by considering the Jacobian matrix J(yk) composed of the derivatives of g evaluated at yk.
g(yk + sk) ≈ g(yk) + J(yk)sk + O(s^2 k)
I neglect the higher order O(s^2 k) terms. Because of this approximation, the sk correction will not make yk + sk a prefect solution. However, the new solution yk + sk will be closer to an ideal solution. Substituting this linear approximation,
Aφk + Ask + g(yk) + J(yk)sk = 0
Solving for sk,
(A + J(φk)) sk = −Aφk − g(yk)
The matrix multiplying sk is known, as are the two quantities on the right. Hence I assume that this system can be solved for sk. Once sk is found, I define yk+1 = yk + sk. yk+1 is then closer to an ideal solution of the differential equation. How- ever, as noted above, yk+1 is not an ideal solution, because of the linear approximation of g(yk + sk). I perform the same process on yk+1 to obtain yk+2, and so on, iterating until the desired accuracy is achieved.
Particular case
We are considering the differential equation,
φ′′^ − sinh φ = 0
The strategy is to break up the interval into a number of small intervals of width h and consider the value of φ at each mesh point. I use the finite differ- ence approximation of the second derivative to write the system of equations explicitly,
φ 0 − 2 φ 1 + φ 2 h^2
− sinh φ 1 = 0 φ 1 − 2 φ 2 + φ 1 h^2
− sinh φ 2 = 0 .. . φn− 2 − 2 φn− 1 + φn h^2
− sinh φ 1 = 0 φn− 1 − 2 φn + φn+ h^2
− sinh φ 1 = 0
Solution
The solutions that I obtain are nearly linear between (− 1 , −2) and (1, 2). I have taken a large number of iterations, and the iterative method appears to be stable.
Code
n = 2^9-1; % Number of subdivisions of the interval h = 2/(n+1); % Size of each subdivision k = 10000; % Number of iterations A = -2eye(n); % Construct the matrix representation of the coefficient matrix for i = 1:n- A(i+1,i) = 1; A(i,i+1) = 1; end ap = -2ones(1, n); % Construct the vector representation of the coefficient matri b = ones(1, n-1); c = ones(1, n-1); phi = linspace(-2+h,2-h,n); % Initial guess for y a = zeros(1, n);% Initialize the other matrices l = zeros(1, n-1); u = zeros(1, n); for t=1:k % Iterative solver, computes LU decomposition l = zeros(1, n-1); % Re-zero the other matrices u = zeros(1, n); y = zeros(1, n); x = zeros(1, n); a = ap-h^2cosh(phi); % A - J(phi_k) f = (-A(phi’))’ + h^2sinh(phi); % -Aphi_k - g(phi_k) f(1) = f(1) + 2; f(n) = f(n) - 2; u(1) = a(1); % Set first value of u for j = 1:n-1 % Compute components of LU decomposition l(j) = c(j)/u(j); u(j+1) = a(j+1) - l(j)*b(j); end y(1) = f(1); % Set first value of y
for j = 2:n % Forwards elimination to find y y(j) = f(j) - l(j-1)y(j-1); end x(n) = y(n)/u(n); % Set last value of x for j = n-1:-1:1 % Backwards elimination to find x x(j) = (y(j) - b(j)x(j+1)) / u(j); end phi = phi + x; % Add iterative correction to phi end close all; % Print and save resulting graphs fp = ’/Users/adrian/Documents/Documents/Academic/School/Classes/Current/Math 128A/ ff = ’-depsc’; fn = ’’; labelstr = ’’; plot([-1:h:1],[-2 phi 2],’b.--’) % Plot the results title([’Iterative solution to D^2\phi - sinh\phi = 0, ’ int2str(n) ’ mesh points, fn = [fp ’Part_1_n_’ int2str(n)]; print(ff,fn);
Part 2
Discussion
Derivation
Multiplying both sides of equation (1) by φ′,
φ′φ′′^ − φ′^ sinh φ = 0
The right side of the above equation is now an anti-derivative of another function. I find this function by integrating. ∫ (^) x
0
φ′(x)φ′′(x) dx −
∫ (^) x
0
φ′(x) sinh φ dx =
∫ (^) x
0
0 dx
I consider each integral separately, ∫ (^) x
0
0 = 0
I make a substitution to perform the first integral,
u = φ′(x) du = φ′′(x)dx
With this substitution, ∫ (^) φ′(x)
φ′(0)
u du =
1 2
u^2 |φ
′(x) φ′(0) =
1 2
(φ′(x))^2 −
1 2
(φ′(0))^2
I also make a substitution to perform the second integral,
u = φ(x) du = φ′(x)dx
With this substitution, ∫ (^) φ(x)
φ(0)
sinh u du = cosh u|φ φ((0)x) = cosh φ(x) − cosh φ(0)
Rearranging terms,
1 2
(φ′(x))^2 − cosh φ(x) =
1 2
(φ′(0))^2 + cosh φ(0)
For any particular value of φ, the quantities on the right are completely determined, and are independent of x. Thus I write them as a constant, E(φ, φ′).
1 2
(φ′(x))^2 − cosh φ(x) = E(φ, φ′)
Solutions
The level curves that I obtain are nearly constant. The variation in the level curves as a function of x decreases as the number of mesh points is increased, indicating that an increased number of points gives a better approximation to the actual function. Note that the MatLab plots below are scaled such that the viewing window is fit to the curve. Although the curves may not look level, notice the small scale on the vertical axis, and that the scale decreases as the number of divisions increases.
Code
Added following the code from part 1,
dif = zeros(1,n); % Finite difference approximation of the first derivative dif(1) = 1./(2h).(phi(2)+2); dif(n) = 1./(2h).(2-phi(n-1)); for j = 2:n- dif(j) = 1./(2h).(phi(j+1)-phi(j-1)); end dplot = 1/2*dif.^2 - cosh(phi); % Level curve to plot figure plot([-1+h:h:1-h],dplot,’g.--’) title([’Approximately level curve, E = 1/2 (D\phi)^2 - cosh\phi , ’ int2str(n) ’ m fn = [fp ’Part_2_n_’ int2str(n)]; print(ff,fn);
Part 3
Discussion
Derivation
Taking E = −1,
−1 =
1 2
(φ′)^2 − cosh φ
⇔
1 2
(φ′)^2 = cosh φ − 1
Now, consider the related hyperbolic function,
sinh^2
φ 2
=
( e
φ 2 − e−^
φ 2
2
) 2
=
( e^2 ·^
φ (^2) − 2 + e^2 ·−^ φ 2
)
4
=
1 2
( eφ^ + e−φ 2
) −
1 2
⇒ 2 sinh^2
φ 2
Using this identity,
1 2
(φ′)^2 = 2 sinh^2
φ 2
⇒ (φ′) 2 = 4 sinh^2 φ ⇒ φ′^ = ±2 sinh φ I solve for φ′^ in the general case to plot level curves, 1 2
(φ′) 2 = E + cosh φ
⇒ φ′^ = ±
√ 2 (E + cosh φ) The phase curves are directed to the left below the φ′^ = 0 axis and to the right above the φ′^ = 0 axis. In the region where φ′^ > 0, φ increases as x is changed, since the derivative is positive. Hence the direction of motion along the curve is to the right, in the direction of increasing φ. Likewise, where φ′^ < 0, changing x causes φ to decrease, so the motion along the phase curve is to the left, in the direction of decreasing φ.
Particular case
Figure 1: Phase space diagram for given differential equation
The phase space diagram for the given differential equation is shown in figure 1. The colored lines correspond to E ≥ −1 and φ′^ > 0. The black lines correspond to solutions with E ≥ −1 and φ′^ < 0, or E < −1. Because of the given boundary conditions, in this case we are primarily interested in the colored curves where E > −1 and φ′^ > 0. The approximate level curves from the iterative solutions are very nearly constant, as noted in part 2 above. The phase space plots of the approximate solutions agree very well with the plotted phase space curves with integer values of E. The approximate solution follows the same shape as the other phase curves, and the approximate phase curve is smooth. It is important that the approximate phase curve does not cross any of the other phase curves, as phase curves can only cross at equilibrium points. The approximate
Results
−2^0 −1.5 −1 −0.5 0 0.5 1 1.5 2
1
2
3
φ
Dφ
Level curve, E = 1/2 (Dφ)^2 − coshφ , 7 mesh points, 10000 iterations
E= E=1 E= E=− E=−1 E(φ , Dφ)
(a)
−2^0 −1.5 −1 −0.5 0 0.5 1 1.5 2
1
2
3
φ
Dφ
Level curve, E = 1/2 (Dφ)^2 − coshφ , 31 mesh points, 10000 iterations
E= E=1 E= E=− E=−1 E(φ , Dφ)
(b)
−2^0 −1.5 −1 −0.5 0 0.5 1 1.5 2
1
2
3
φ
Dφ
Level curve, E = 1/2 (Dφ)^2 − coshφ , 127 mesh points, 10000 iterations
E= E= E=0 E=− E=− E(φ , Dφ)
(c)
−2^0 −1.5 −1 −0.5 0 0.5 1 1.5 2
1
2
3
φ
Dφ
Level curve, E = 1/2 (Dφ)^2 − coshφ , 511 mesh points, 10000 iterations
E= E= E=0 E=− E=− E(φ , Dφ)
(d)
