source("http://macosa.dima.unige.it/r.R") # If I have not already loaded the library ---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------S09Differential equations (1st and 2nd order). Direction fields# Graphic solution of1st order differential equationsby the Euler method:#soledif(x0,y0,xf,N,col) traces with N steps the polygonal approximation of the # solution of y' = f(x,y) (that passes through (x0,y0) and goes from x0 to xf); # f is in Dy. #soledif1(x0,y0,xf,N,col) traces itslim. #diredif(a,b,c,d, m,n) draws (for a<=x<=b, c<=y<=d) thedirection fieldwith m rows # and n columns # [you have to useyandDyto describe the differential equation; # instead of x you can use another variable]# It is always advisable to trace the directional field to deal with problems relating # to differential equations. It is easy to make mistakes in solving differential # equations and mathematical software (like Wolfram Alpha) sometimes gives incorrect # solutions.#Example:y' = x-ysuch that (x0,y0) = (-1,3), for -1 <= x <= 5Dy= function(x,y) x-yPlane(-1,5, -2,5)soledif(-1,3, 5, 10, "brown") soledif(-1,3, 5, 1e3, "blue")# I understand that with N = 1000 the graph has stabilized. Typically, N = 1e4soledif(4,2.5, -1, 1e4, "red")# It's enough; if the graph has a lot of oscillation it is necessary to increase N # If I want to see the directional field for 20 rows and 25 columns I dodiredif(-1,5, -2,5, 25,20)# If I just want to see this I do not put the "soledif" lines.# Other examples for the same equation (right figure):Plane(-1,5, -2,5) soledif(1,0, 5, 1e4, "brown"); soledif(1,0, -1, 1e4, "blue") soledif(2,1.5, 5, 1e4, "brown"); soledif(2,1.5, -1, 1e4, "blue") soledif(3,1, 5, 1e4, "brown"); soledif(3,1, -1, 1e4, "blue")#soledifv(x0,y0, p,N) gives thevalue of y(p)approximated with N steps P=soledifv(3,1, 4,1e5); P; P=soledifv(3,1, 4,1e6); P#2.632122 2.632121 POINT(p,P, "red") # Theredpoint in the picture#soledifp(x0,y0,xf,N,col) tracks only N + 1 points without linking themPlane(-1,5, -2,5) soledifp(-1,3, 5, 10, "black") Plane(-1,5, -2,5) soledif(-1,3, 5, 10, "red") soledifp(-1,3, 5, 10, "black")# In this case (y'=x-y) I know how to solve the equation; the solution is: f = function(x) x-1+(C+1)*exp(-x) # # The solution such that f(3)=1: # f(3) = 1 implies 2+(C+1)*exp(-3) = 1, C = -1/exp(-3)-1 C = -1/exp(-3)-1 # I can find the value found before: f(4)#2.632121 OK # If I want I can compare the graphs of the two solutions Plane(-1,5, -2,5) POINT(3,1,"black"); soledif(3,1, 5, 1e4, "brown") graph1(f, -1,5, "blue")# I can find the solution withWolframAlpha: # y'(x) = x-y(x), y(3) = 1, y(4) gives 3-1/exp(1) (= 2.632121)#Another example:y'(x) = 3·yDy = function(x,y) 3*(y^2)^(1/3) Plane(-2,2, -3,3) diredif(-2,2, -3,3, 20,20) soledif(1,1, 2, 1e4, "seagreen") soledif(1,1, -2, 1e4, "seagreen") soledif(0,0, 2, 1e4, "blue") soledif(0,0, -2, 1e4, "brown") soledif(-2,-1, -1, 1e4, "brown")^{2/3}(x)# It is easy verify that y(x) = x^{3}is a solution such thaty(1)=1, but this is not the # only solution! (gives only this solution!). # This example shows how useful it is to trace the directional field.WolframAlpha#A third example:y'(x) = x/y(x)such that y(-3)=2; we want find the value of y(-1.5). Dy = function(x,y) x/y Plane(-3,1, -2,2); diredif(-3,1, -2,2, 22,22) POINT(-3,2, "black") #From the directional field(picture on the left) I can deduce that there is no solution. # Let's face the problem. soledif(-3,2, -1.5, 1e5, "seagreen") # I obtain the green graph (picture on the right); and if I look for the value of the # solution I find: soledifv(-3,2, -1.5, 1e5)#-2.525611 ?!? # I can usesoledifvcthat prints the coordinatevalues of the first four couples of # points where theslopeof the approximate curve-solution:changes the signsoledifvc(-3,2, -1.5, 1e6) # -2.236062 -0.00046688 -2.23606 0.006717179 # -2.23606 0.006717179 -2.236059 0.006217848 # -2.236048 -0.0006352839 -2.236047 0.00464436 # -2.236047 0.00464436 -2.236046 0.003922179 # I can improve the approximation of the value: soledifvc(-3,2, -2.23, 1e6)#-2.236065 -0.001458657 -2.236064 -0.0002782769 POINT(-2.236,0, "brown") # Another question. If there are somepoints where the equation is non defined, it may # happen that they are not identified by the program. # An example with the same equation: PLANE(-3,3, -3,3) diredif(-3,3, -3,3, 20,20) soledif(-3,3, 3, 1e4, "brown") soledif(-3,-3, 3, 1001, "blue") # It may happen that the program overrides x=0 where the equation is not defined.#A fourth example:y'(x) = (x+y(x))/(x-y(x))such that y(-2)=1. Dy = function(x,y) (x+y)/(x-y) BF=2.7; HF=2.7; Plane(-3,3, -3,3) diredif(-3,3, -3,3, 20,20) soledif(-2,1, -3, 1e5, "red") soledif(-2,1, 3, 1e5, "red") BF=2.7; HF=2.7; Plane(-3,3, -3,3) soledif(-2,1, -3, 1e5, "red") soledifvc(-2,1, 3, 1e6)#-1.146135 1.146135 -1.14613 1.146135#0.238275 0.2370324 0.23828 0.238945 # 0.23828 0.238945 0.238285 0.2353571 # 0.238285 0.2353571 0.23829 0.236166 # The first point where the slope changes the sign is amaximum pointsoledif(-2,1, 0.2384, 1e5, "red") POINT(-1.146,1.146, "brown"); POINT(0.2384,0.2345, "blue")#A fifth example: integral curves ofy'(x) + x*y(x)^2/(1+x^2*y(x)) = 0Dy = function(x,y) -x*y^2/(1+x^2*y) BF=5; HF=4; Plane(-2.5,2.5, -4,3) diredif(-2.5,2.5, -4,3, 21,21) # through the point (0,2) soledif(0,2, -3, 1e5, "brown"); soledif(0,2, 3, 1e5, "brown") soledif(0,1, 3, 1e5, "brown"); soledif(0,1, -3, 1e5, "brown") # through the point (0,-1) soledifvc(0,-1, -3, 1e6) # 0 -1 -3e-06 -1 # -0.707118 -2.00314 -0.707121 -1.997825 # -0.707121 -1.997825 -0.707124 -2.005909 # -0.707124 -2.005909 -0.707127 -2.003067 soledif(0,-1, -0.70711, 1e5, "seagreen"); soledif(0,-1, 0.70711, 1e5, "seagreen") soledif(-0.7073, -2, -0.1, 1e5, "seagreen"); soledif(0.7073, -2, 0.1, 1e5, "seagreen") # through the point (0,-1/2) soledifvc(0,-0.5, -3, 1e6) # 0 -0.5 -3e-06 -0.5 # -1.000011 -1.000668 -1.000014 -0.9963146 # -1.000014 -0.9963146 -1.000017 -0.9971288 # -1.000026 -1.031578 -1.000029 -1.031477 soledif(0,-1/2, -1, 1e5, "seagreen"); soledif(0,-1/2, 1, 1e5, "seagreen") soledif(-1.0001, -1.0001, -0.5, 1e5, "seagreen"); soledif(1.0001, -1.0001, 0.5, 1e5, "seagreen") # through the points (-2,-1), (2,-1) soledif(-2,-1, -3, 1e5, "blue"); soledif(-2,-1, -0.1, 1e5, "blue") soledif(2,-1, 3, 1e5, "blue"); soledif(2,-1, 0.1, 1e5, "blue") # Of course, the simplest differential equations areindefinite integrals, which can be # studied with simpler methods (see). To trace the directional field it is useful to # solve other simple differential equations, ot to better understand the solution. An # example. The differential model f'(x) = f(x) [that is the differential equationy' =#y] is easy to solve: I know that y=exp(x) is a solution, and that [because D_{x}(k*f(x)) # = k*D_{x}f(x)] y(x) = k*exp(x) is also a solution (for every number k). Let's check. Dy = function(x,y) y Plane(-2,2, -5,5); diredif(-2,2, -5,5, 15,15) soledif(0,1, 2, 1e5, "brown"); soledif(0,1, -2, 1e5, "brown") soledif(-0.5,-2, 2, 1e5, "seagreen"); soledif(-0.5,-2, -2, 1e5, "seagreen") # Solving the differential equation: # k*exp(1)=-2 => k=-2/exp(1) g=function(x) -2/exp(1)*exp(x); POINT(1,-2,"brown") graph1(g, -2,2, "black") # Aproblem: # A run-down condenser has capitance C = 80 μF and is connected to a battery of 45 V # with a resistor of resistance R = 500 Ω. What is the final charge on the condenser?V=45; C=80*10^-6; R=500 # y is Q (= C*V); Q(0)=0 is the initial condition Dy = function(t,y) (V-y/C)/R# Try looking for the solution with a particular graphic scale.Plane(0,10, 0,1); soledif1(0,0, 10, 1e5, "blue") # Then, based on the results, I modify it. Plane(0,10, 0,0.05); soledif1(0,0, 10, 1e5, "blue")# I modify it.Plane(0,0.35, 0,0.004); soledif1(0,0, 0.35, 1e5, "blue")# I look for the limit (from 0 to 0.35 or - better - from 0 to 1) soledifv(0,0, 1, 1e5)# 0.0036coldash="brown"; segm(-1,0.0036, 0.4, 0.0036, 0) abovex("t"); abovey("Q")# It is not difficult to solve this equation "by hand". # Let's see how to solve it withWolfram Alpha# I put Q(t)/C + R*dQ(t)/dt = V, Q(0)=0 # [or dy(t)/dt = (V-y(t)/C)/R, y(0)=0 ]# I obtain:Q(t) = C*V*(1-exp(-t/(C*R)))# Verification:F = function(t) C*V*( 1-exp(-t/(C*R)) ); graph1(F, 0,0.35, "brown")# I get a graph that overlaps with the previous one! # We can verify the solutions ofpartial differential equations. One example. f = function(x,y) log(x^2+y^2) # It is a solution ofderiv2(f,"x"); deriv2(f,"y") # 2/(x^2 + y^2) - 2 * x * (2 * x)/(x^2 + y^2)^2 A # 2/(x^2 + y^2) - 2 * y * (2 * y)/(x^2 + y^2)^2 B # It is easy to check algebraically that A+B = 0. Let's check it numerically: dxxf = function(x,y) eval(deriv2(f,"x")); dyyf = function(x,y) eval(deriv2(f,"y")) g = function(x,y) dxxf(x,y)+dyyf(x,y) x=runif(4, 0,20); y=runif(4, 0,20); g(x,y) # -4.336809e-19 0.000000e+00 -8.673617e-19 1.734723e-18 OK: practically 0 # Other solutions, for any k: f = function(x,y) exp(k*x)*cos(k*y) # I take for example: k = 3/7 deriv2(f,"x"); deriv2(f,"y") # exp(k * x) * k * k * cos(k * y) A # -(exp(k * x) * (cos(k * y) * k * k)) B # Numerical verification (instead of the algebraic one) that A + B = 0: dxxf = function(x,y) eval(deriv2(f,"x")); dyyf = function(x,y) eval(deriv2(f,"y")) g = function(x,y) dxxf(x,y)+dyyf(x,y) x=runif(4, 0,20); y=runif(4, 0,20); g(x,y) # -2.842171e-14 5.551115e-17 2.842171e-14 0.000000e+00 OK: practically 0 # For∂/∂x ∂/∂x f(x,y) + ∂/∂y ∂/∂y f(x,y) = 0 seesecond-orderdifferential equationshere.Other examples of use