# Warning: These snips of code made for y'=1-x-y, y(0)=3. # Code computes approx values for y(0.1) to y(1). # 'Dots' is the list of dots for connect-the-dots graphics. # ======================================== # Euler. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=200: # Group 2, repeat n times. Euler's method for i from 1 to n do Y:=y0+h*f(x0,y0); x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display relevant dots and plot. Exact:=x->2-x+exp(-x); P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)): m:=40:X:=[seq(1+m*j,j=0..n/m)]: # List of relevant indices print("Dots"),seq(Dots[k],k=X); print("Exact"),seq(Exact(Dots[k][1]),k=X); print("Error"),seq(P(Exact(Dots[k][1]),Dots[k][2]),k=X); plot([Dots]); # ======================================== # Heun. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=200: # Group 2, repeat n times. Heun method. for i from 1 to n do Y1:=y0+h*f(x0,y0); Y:=y0+h*(f(x0,y0)+f(x0+h,Y1))/2: x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display relevant dots and plot. Dots[1],Dots[2],seq(Dots[1+40*j],j=1..n/40); plot([Dots]); # ======================================== # RK4. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=100: # Group 2, repeat n times. RK4 method. for i from 1 to n do k1:=h*f(x0,y0): k2:=h*f(x0+h/2,y0+k1/2): k3:=h*f(x0+h/2,y0+k2/2): k4:=h*f(x0+h,y0+k3): Y:=y0+(k1+2*k2+2*k3+k4)/6: x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display some dots and plot. Dots[1],Dots[2],Dots[101]; plot([Dots]); # Code snips for exact/error reports # ========================================= # Making multiple curves on one plot # ======================================== Exact:=x->2-x+exp(-x); # An exact solution plot({Exact(x),[Dots]},x=0..1/2); # plot exact and approx solutions # ======================================== # How to create a Dots table for the exact solution # ======================================== Exact:= x -> 2-x+exp(-x):n:=10: ExactDots:=seq([Dots[j][1],Exact(Dots[j][1])],j=1..n+1); # ======================================== # How to define and print percentage relative error: # ======================================== P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)); ExactVal:=Exact(Dots[11][1]): # Compute exact y-value for x=0.5 ApproxVal:=Dots[11][2]: # Get Euler approx y-value for x=0.5 P(ExactVal),ApproxVal); # print percent relative error # ======================================== # How to create a Dots table for percentage error # ======================================== P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)); Pdots:=seq([Dots[j][1],P(Exact(Dots[j][1]),Dots[j][2])],j=1..11); # ========================================= # Printing results and tables # Make tables with a pencil, it saves time. # ======================================== # To extract and print items 1,101,201,1001 from a list: Dots1:=Dots[1],Dots[101],Dots[201],Dots[1001];