• API
• FAQ
• Tools
• Archive
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.

Top