Error converting function handle to double?

Hey guys, any idea why this won't run? I attached the code below as well as the original problem. Thanks!

% HW 9

% define the mesh in space
dx = 1e-2;
x = 0:dx:10;
x = x';

% define the mesh in time
dt = 1e-4;
t = 0:dt:2;

c = 1/2;
% display('this number should be between -1 and 1 for stability:')
mu = c*dt/dx;

% choose the wave number of the initial data and give its decay rate
u = zeros(N+1,M+1);
u(:,1) = @(t) exp(-20*((x-3)^2));
u_exact(:,1) = u(:,1);

% I want to do the Lax Wendroff scheme:
%
% u_new(j) = u_old(j) - (mu/2)*(u_old(j+1)-u_old(j-1)) 
%                           + (mu^2/2)*(u_old(j-1)-2*u_old(j)+u_old(j+1))

for j=1:M
    for k=2:N
        u(k,j+1) = u(k,j) - (mu/2)*(u(k+1,j)-u(k-1,j)) + (mu^2/2)*(u(k+1,j)-2*u(k,j)+u(k-1,j));
    end
    % I code in the exact values at the endpoints.
    u(1,j+1)=1;
    u(N+1,j+1)=0;
    X = x-c*t(j+1);
    u_exact(:,j+1) = @(t) 0.4*sin(2*t);
end

Answer 1

[Geoff Hayes]

Pizza Pie - the full error message is

 The following error occurred converting from function_handle to double:
 Error using double
 Conversion to double from function_handle is not possible.

 Error in ***** (line 18)
 u(:,1) = @(t) exp(-20*((x-3)^2));

You have defined u as a matrix of zeros (doubles) with the line of code

 u = zeros(N+1,M+1);

and then you are trying to assign a anonymous function handle to all elements in the first column of u. Hence the error.

I suspect that you want to evaluate the expression instead though it is unclear how t should be used in it. If, for the moment, we ignore the t then something like the following will work

 u(:,1) = exp(-20*((x-3).^2));

(The above assumes that N is 1000.)

What is the expression that you wish to evaluate? How should t be used?

[Pizza Pie]

I changed it to the following and am still getting an error

% define the mesh in space
dx = 1e-2;
x = 0:dx:10;
x = x';

% define the mesh in time
dt = 1e-4;
t = 0:dt:2;

c = 1/2;
% display('this number should be between -1 and 1 for stability:')
mu = c*dt/dx;

% choose the wave number of the initial data and give its decay rate
u = zeros(N+1,M+1);
u(:,1) = exp(-20*((x-3).^2));
u_exact(:,1) = u(:,1);

% I want to do the Lax Wendroff scheme:
%
% u_new(j) = u_old(j) - (mu/2)*(u_old(j+1)-u_old(j-1)) 
%                           + (mu^2/2)*(u_old(j-1)-2*u_old(j)+u_old(j+1))

for j=1:M
    for k=2:N
        u(k,j+1) = u(k,j) - (mu/2)*(u(k+1,j)-u(k-1,j)) + (mu^2/2)*(u(k+1,j)-2*u(k,j)+u(k-1,j));
    end
    % I code in the exact values at the endpoints.
    u(1,j+1)=1;
    u(N+1,j+1)=0;
    X = x-c*t(j+1);
    u_exact(:,j+1) = @(t) 0.4*sin(2*t);
end

I attached the full problem image above if you wanted to see the goal of the problem. Maybe I am doing something else wrong?

[Geoff Hayes]

There is no image attached.

Please copy and paste the full error message. Also, what are N and M?

[Pizza Pie]

Image has been attached now. Apologies

Answer 2

[Star Strider]

The problem is with your anonymous functions. I created versions that work, and made what seem to me to be the appropriate substitutions in your ‘u’ and ‘u_exact’ assignments. (You didn’t supply ‘N’ and ‘M’, so I made up some values for them to test the code.) I will leave you to determine if your code with my changes produces the correct results:

N = 3;
M = 4;
% define the mesh in space
dx = 1e-2;
x = 0:dx:10;
x = x';

% define the mesh in time
dt = 1e-4;
t = 0:dt:2;

c = 1/2;
% display('this number should be between -1 and 1 for stability:')
mu = c*dt/dx;

% choose the wave number of the initial data and give its decay rate
u = zeros(N+1,M+1);
u_fcn = @(x) exp(-20*((x-3).^2));                                       % <— CONVERT TO ANONYMOUS FUNCTION
u_exact_fcn = @(t) 0.4*sin(2*t);                                        % <— CONVERT TO ANONYMOUS FUNCTION
u_exact(:,1) = u_fcn(u(:,1));

% I want to do the Lax Wendroff scheme:
%
% u_new(j) = u_old(j) - (mu/2)*(u_old(j+1)-u_old(j-1)) 
%                           + (mu^2/2)*(u_old(j-1)-2*u_old(j)+u_old(j+1))

for j=1:M
    for k=2:N
        u(k,j+1) = u_fcn(u(k,j)) - (mu/2)*(u_fcn(u(k+1,j))-u_fcn(u(k-1,j))) + (mu^2/2)*(u_fcn(u(k+1,j))-2*u_fcn(u(k,j))+u_fcn(u(k-1,j)));
    end
    % I code in the exact values at the endpoints.
    u(1,j+1)=1;
    u(N+1,j+1)=0;
    X = x-c*t(j+1);
    u_exact(:,j+1) = u_exact_fcn(t(j));
end

[Pizza Pie]

This isn't giving me any error messages so far, however - how would i go about plotting this every 10% of the way? I can't use the normal plot function like this so I'm guessing it has to be in the for loop somewhere?

[Star Strider]

What do you want to plot?

Also, does the notation in the screen shot image indicate that this is a homegeneous partial differntial equation?