Wave-equation/wave1d.m at main · SarahShakim/Wave-equation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
% Solving the wave equation
%
% Parameters
% ==========
%    c = speed of the wave
%    x_int = x interval
%    n_x = number of values to split x interval into
%    t_int = t interval
%    n_t = number of values to split t interval into
%    u_init =  initial value function
%    du_init = intital rate of change of function.
%    u_bndry = boundary condition function
%
% Return Values
% =============
%    x_out = The vector x is the boundary of values
%    t_out = The vector t is the boundaries for the matrix
%    U_out =  matrix giving ghe temperature values for each point in space
%    in a certain time

function [x_out, t_out, U_out] = wave1d( c, x_int,  n_x, t_int, n_t, u_init, du_init, u_bndry )

    % ARGUMENT CHECKING


    %initializing the variables for the boundaries:
    a = x_int(1);
    b = x_int(2);

    %initializing the boundaries for
    t_initial = t_int(1);
    t_final = t_int(2);


    dt = (t_final - t_initial)/(n_t - 1);

    %determining the h value
    h = (b - a)/(n_x - 1);

    %determining the value of r
    r = (c*dt/h)^2;

      %initializing the matrix of zeroes
    U_mat = zeros(n_x, n_t);


    % ERROR CHECKING
    % Calculate the error ratio (c * del_t / h)^2 to see if it is greater or equal to 1
        err = (c*dt/h)^2;

        if err>= 1

            noft = ceil((c*(t_final-t_initial)/h) + 1);

            throw( MException( 'MATLAB:invalid_argument', ...
                'The ratio (c * del_t / h)^2 = %f >= 1, consider using n_t = %d', ...
                err, noft ) );
        end

    %creating a t vector using n_t
    t_vec = linspace(t_initial, t_final, n_t);

    %creating an x vector using n_x
    x_1 = linspace(a, b, n_x)';
    x_vec = u_init(x_1);

    % Assign intial values using x vector
    U_mat(:,1) = x_vec;

    % putting the boundaries into the vector
    bound = u_bndry(t_vec);
    top_bound = bound(1);
    bot_bound = bound(2);


    % Assign boundary space values over time
    U_mat(1, 2:end) = top_bound;
    U_mat(end, 2:end) = bot_bound;

    % SOLVING
    % used the finite difference fomula to find every value in the matrix
    for t = 2:n_t

        % Special case for t = t_2
        % Use Euler's approximation.
        if t == 2
            u_x2 = U_mat(2:end - 1, 1) + dt*(du_init(x_vec(2:end - 1)));
            U_mat(2:end - 1, 2) = u_x2;


            if isnan(U_mat(1, 2))
                % Solve using forward and backward divided difference approximation
                U_mat(1, 2) = (4/3)*U_mat(2, 2) - (1/3)*U_mat(3, 2);
            end
            if isnan(U_mat(end, 2))
                % Solve using forward and backward divided difference approximation
                U_mat(end, 2) = (4/3)*U_mat(end - 1, 2) - (1/3)*U_mat(end - 2, 2);
            end


        else
            u_i_k = U_mat(:, t-1);
            u_i_k1 = 2*(U_mat(2:end - 1, t-1)) - (U_mat(2:end - 1, t-2)) + r.*(diff(u_i_k, 2));

            U_mat(2:end - 1, t) = u_i_k1;


            if isnan(U_mat(1, t))
                %using forward and backward divided difference approximation
                U_mat(1, t) = (4/3)*U_mat(2, t) - (1/3)*U_mat(3, t);
            end
            if isnan(U_mat(end, t))
                %using forward and backward divided difference approximation
                U_mat(end, t) = (4/3)*U_mat(end - 1, t) - (1/3)*U_mat(end - 2, t);
            end
        end
    end

    x_out = x_1;
    t_out = t_vec;
    U_out = U_mat;

end