SHARE
TWEET

Untitled

a guest Jun 7th, 2019 93 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. t0 = 0;
  2. x0 = 1;
  3. v0 = 0;
  4. T = 40;
  5. A = 0;
  6. B = 0;
  7. C = 0.3;
  8. a = 1;
  9. b = 1;
  10. c = 0;
  11. N = 10000;
  12.  
  13. epsilon = 1e-4;
  14.  
  15. f = @(t, x) [(a - c * x(2) ^ 2) * x(1) - b * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];
  16. f_b_plus = @(t, x) [(a - c * x(2) ^ 2) * x(1) - (b + epsilon) * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];
  17. f_b_minus = @(t, x) [(a - c * x(2) ^ 2) * x(1) - (b - epsilon) * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];
  18.  
  19. s_value = sqrt(a*a - 4 * b);
  20.  
  21. x_correct_f = @(t) exp(t/2 * (a - s_value)) * (exp(t* s_value) * (s_value - a) + s_value + a) / (2 * s_value);
  22. v_correct_f = @(t) -b * exp(t/2 * (a - s_value)) * ((exp(t* s_value) - 1));
  23.  
  24. t = linspace(t0, T, N + 1);
  25. xvs = zeros(2, N + 1);
  26. vs = zeros(1, N + 1);
  27. xs = zeros(1, N + 1);
  28. xvs(1,1) = v0;
  29. xvs(2,1) = x0;
  30.  
  31. for i = 1:N
  32.   xvs(:,i + 1) = schemat(f, t(i), xvs(:,i), t(i + 1), 1);
  33.   vs(i + 1) = v_correct_f(t(i));
  34.   xs(i + 1) = x_correct_f(t(i));
  35. endfor
  36.  
  37.  
  38. mean_abs_diff_x = max(abs(xvs(2, :) - xs))
  39. mean_abs_diff_v = max(abs(xvs(1, :) - vs))
  40. printf("mean abs diff with x %f,\nmean abs diff with v %f\n", mean_abs_diff_x, mean_abs_diff_v)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top