Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Crank-Nicolson method

Subject: Crank-Nicolson method

From: Scott

Date: 22 Aug, 2010 12:21:04

Message: 1 of 17

I'm trying to follow an example in a MATLab textbook. It has the following code which I have simply repeated. After the code it says: "the following MATLab function heat_crank.m finds the solution of the heat equation using the Crank-Nicolson method. Inputs are function f; boundary condition c1, c2; endpoint L; maximum time T; step sizes h and k; constant alpha. The input function f(x) should be defined as an m-file."

I have values:
c1=2;
c2=2;
L=2;
T=0.5;
h=0.2;
k=0.02;
alpha=1;
and f(x)=2*(x-1)^3

What does it mean the input function f(x) should be defined as an m-file??
Do I need to create an m-file with the above values in it and then say somewhere inside heat_crank.m to look at this other m-file for the values? Or do I just put these values into the heat_crank.m file somewhere at the start?

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 22 Aug, 2010 15:23:42

Message: 2 of 17

On 22/08/10 7:21 AM, Scott wrote:
> I'm trying to follow an example in a MATLab textbook. It has the
> following code which I have simply repeated. After the code it says:
> "the following MATLab function heat_crank.m finds the solution of the
> heat equation using the Crank-Nicolson method. Inputs are function f;
> boundary condition c1, c2; endpoint L; maximum time T; step sizes h and
> k; constant alpha. The input function f(x) should be defined as an m-file."

> What does it mean the input function f(x) should be defined as an m-file??
> Do I need to create an m-file with the above values in it and then say
> somewhere inside heat_crank.m to look at this other m-file for the
> values? Or do I just put these values into the heat_crank.m file
> somewhere at the start?

It means you should create a file named f.m that contains

function y = f(x)
   y = 2*(x-1)^3;
end

As long as that file is on the regular matlab search path, you will not
need to tell the heat_crank.m program where to look for f.m . However,
depending exactly what kind of input heat_crank is looking for, in the
place where it is expecting you to pass f, you may have to pass 'f' (the
string containing the letter f), or you may have to pass @f (ampersand
followed by f, with no quotation marks; this will also depend upon your
MATLAB version.


The documentation for heat_crank.m is... incomplete... when it tells you
to do things this way, but it is Sunday morning here and you didn't need
to repost your question so early, so I'm just giving you the answer in
exact accordance with the documentation. You will either eventually
learn better or you won't.

Subject: Crank-Nicolson method

From: Scott

Date: 22 Aug, 2010 15:37:06

Message: 3 of 17

Ok, so I tried that but it hasn't worked because I haven't defined values for c1,c2 etc etc. In which file do I do this??

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 22 Aug, 2010 16:29:42

Message: 4 of 17

On 22/08/10 10:37 AM, Scott wrote:
> Ok, so I tried that but it hasn't worked because I haven't defined
> values for c1,c2 etc etc. In which file do I do this??

You need to spend some time reading the Getting Started documentation.

Subject: Crank-Nicolson method

From: ImageAnalyst

Date: 22 Aug, 2010 17:11:48

Message: 5 of 17

On Aug 22, 11:37 am, "Scott " <boofheadf...@hotmail.com> wrote:
> Ok, so I tried that but it hasn't worked because I haven't defined values for c1,c2 etc etc. In which file do I do this??

You do that in your main function. For example:
In test_f.m:

function test_f
    c1=2;
    c2=2;
    L=2;
    T=0.5;
    h=0.2;
    k=0.02;
    alpha=1;
    y = f(1:20);
    plot(y);

And in f.m:
   function y = f(x)
        y=2*(x-1).^3; % Note .^ not just ^

Subject: Crank-Nicolson method

From: Scott

Date: 23 Aug, 2010 01:00:25

Message: 6 of 17

Hi, instead of putting the function for f in a separate m-file I tried putting it in the main file. I also tried to define c1,c2,L etc values in the main function. My script is below. When I try to run the script it says this:

??? Input argument "f" is undefined.
Error in ==> heat_crank at 16
    u(i)=feval(f,(i-1)*h);
Error in ==> run at 57
          evalin('caller', [s ';']);

This is my script:

function heat_crank(f,c1,c2,L,T,h,k,alpha)
c1=2;
c2=2;
L=2;
T=0.5;
h=0.02;
k=0.02;
alpha=1;

n=L/h;
m=T/k;
lambda=alpha^2*k/(h^2);
z=0:h:L;

for i=1:n+1
    u(i)=feval(f,(i-1)*h);
end

for i=2:n
    if(i~=n)
        a(i)=-lambda;
    end
    b(i)=2*(1+lambda);
    if (i~=n)
        c(i)=-lambda;
    end
end

bb=b;
for j=1:m
    t=j*k;
    
    for i=2:n
        d(i)=lambda*u(i-1)+2*(1-lambda)*u(i)+lambda*u(i+1);
    end
    y(n+1)=c1;
    y(1)=c2;
    for i=3:n
        ymult=a(i-1)/bb(i-1);
        bb(i)=bb(i)-ymult*c(i-1);
        d(i)=d(i)-ymult*d(i-1);
    end
    y(n)=d(n)/bb(n);
    for i=n-1:-1:2
        y(i)=(d(i)-c(i)*y(i+1))/bb(i);
    end
    u=y;
    bb=b;
end
        
function y = f(x)
    y = 2*(x-1).^3;
    return

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 23 Aug, 2010 02:09:53

Message: 7 of 17

On 22/08/10 8:00 PM, Scott wrote:
> Hi, instead of putting the function for f in a separate m-file I tried
> putting it in the main file. I also tried to define c1,c2,L etc values
> in the main function. My script is below. When I try to run the script
> it says this:
> ??? Input argument "f" is undefined.
> Error in ==> heat_crank at 16

Either you did not read the Getting Started section or else you did not
understand it.

Subject: Crank-Nicolson method

From: Scott

Date: 23 Aug, 2010 02:24:05

Message: 8 of 17

Clearly didn't understand it,

I presumed I may have used the y variable in too many instances so I relabled it u because that's what the function f(x) is, but still no luck. When I try to end the function there is also a problem there because you can't end that function for some reason... for sure I know it will be as easy as changing one or two lines but I'm stuck as to which.

Scott

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 23 Aug, 2010 02:38:37

Message: 9 of 17

On 22/08/10 9:24 PM, Scott wrote:
> Clearly didn't understand it,
>
> I presumed I may have used the y variable in too many instances so I
> relabled it u because that's what the function f(x) is, but still no
> luck. When I try to end the function there is also a problem there
> because you can't end that function for some reason... for sure I know
> it will be as easy as changing one or two lines but I'm stuck as to which.

function A_Blatant_Hint
   a = 1;
   b = 2;
   c = cranky(@f, a, b)
end

Subject: Crank-Nicolson method

From: Scott

Date: 23 Aug, 2010 02:58:05

Message: 10 of 17

> function A_Blatant_Hint
> a = 1;
> b = 2;
> c = cranky(@f, a, b)
> end

Walter, I'm sorry this is so obvious but it's my first time using MATLab. As easy as it is to you I didn't understand that last message one bit. I have:

function u = f(x)
    u = 2*(x-1).^3;
end

and I don't understand what is wrong with that. u is my function which is 2(x-1)^3. I wrote it as above. I assigned values to c1,c2, up to alpha at the beginning on the script. now, down the bottom of my script after the function heat_crank I have created a new function which I want to be f, that will describe the function f(x) to equal 2(x-1)^3. It doesn't work, there is always an error. The person previously told me to write what I did and I put it inside the same m-file but still nothing...

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 23 Aug, 2010 03:17:56

Message: 11 of 17

On 22/08/10 9:58 PM, Scott wrote:
> The person previously told me to
> write what I did

No they didn't: they said to put things in your main file. heat_crank is
not your main file.

 > and I put it inside the same m-file but still nothing...

Let's start with something simpler

function r = Another_blatant_Hint(p, q)
   r = p + q;
end

Now, if you were to store that in the file Another_blatant_hint.m
then how would you use it to add (say) 19 and 62 ? Which buttons would
you press, or what would you type at the command line?

Subject: Crank-Nicolson method

From: Marc

Date: 23 Aug, 2010 03:50:08

Message: 12 of 17

Scott

Got your script to work, should be pretty straight forward although I am not sure if the math is right or wrong.

Since you are defining c1, c2, etc. etc. you did not need to have that in your leading line. In your original code you had left out an end for the 1st function line heat_crank.

In the below code, I commented out the function heat_crank111 call along with the "sub" function food and their corresponding code so that the various values would end up in my workspace just to make sure it did something....

Not sure if there was anything wrong with your use of feval, I just do not typically use feval for something that straight forward so I went to my comfort zone.

Hope this helps.

%function heat_crank111
c1=2;
c2=2;
L=2;
T=0.5;
h=0.02;
k=0.02;
alpha=1;

n=L/h;
m=T/k;
lambda=alpha^2*k/(h^2);
z=0:h:L;

for i=1:n+1
    u(i)=food((i-1)*h);
end

for i=2:n
    if(i~=n)
        a(i)=-lambda;
    end
    b(i)=2*(1+lambda);
    if (i~=n)
        c(i)=-lambda;
    end
end

bb=b;
for j=1:m
    t=j*k;
    
    for i=2:n
        d(i)=lambda*u(i-1)+2*(1-lambda)*u(i)+lambda*u(i+1);
    end
    y(n+1)=c1;
    y(1)=c2;
    for i=3:n
        ymult=a(i-1)/bb(i-1);
        bb(i)=bb(i)-ymult*c(i-1);
        d(i)=d(i)-ymult*d(i-1);
    end
    y(n)=d(n)/bb(n);
    for i=n-1:-1:2
        y(i)=(d(i)-c(i)*y(i+1))/bb(i);
    end
    u=y;
    bb=b;
end

%end
        
%function f = food(x)
% f = 2*(x-1).^3;
%end

-------------------------------
Save this as a seperate m-file, food.m or whatever you would like to call your function
-----------------------------
function f = food(x)
    f = 2*(x-1).^3;
end

Subject: Crank-Nicolson method

From: Scott

Date: 23 Aug, 2010 04:52:03

Message: 13 of 17

I took away the header function heat_crank and made a new m-file called f_x with the function in it and ran heat_crank and it appears to have worked. I have to show how the temperature (z-axis) appears over time and length so I want to make a mesh grid of the values of temperature over the time 0<t<0.5 and length 0<x<2. I can't get this to work. I'm not sure if all the temperature values have been calculated or if it's just calculated temperature for a certain x or t value...

Any suggestions?

Subject: Crank-Nicolson method

From: Walter Roberson

Date: 23 Aug, 2010 05:39:55

Message: 14 of 17

On 22/08/10 11:52 PM, Scott wrote:
> I took away the header function heat_crank and made a new m-file called
> f_x with the function in it and ran heat_crank and it appears to have
> worked. I have to show how the temperature (z-axis) appears over time
> and length so I want to make a mesh grid of the values of temperature
> over the time 0<t<0.5 and length 0<x<2. I can't get this to work. I'm
> not sure if all the temperature values have been calculated or if it's
> just calculated temperature for a certain x or t value...
>
> Any suggestions?

My suggestion is to look at the size() of the output first, to verify
that one of the dimensions is the number of steps you allowed for the
temperature, and that the other dimension is the number of steps you
allowed for the length. Once you have verified that your program is
calculation the correct number of outputs, use surf() to plot the
output: if all of the outputs in one direction are exactly the same as
the others, then you need to double-check your calculation. (Sometimes
all of the outputs being the same is correct but unexpected -- sometimes
it is just the way the math works out. But _usually_ it indicates an
error in the calculation method.)

Subject: Crank-Nicolson method

From: Scott

Date: 23 Aug, 2010 08:27:04

Message: 15 of 17

I've done a complete rewrite and I have this code. It runs fine except it is not quite calculating the correct values. I think there is something wrong with my code for uu or somewhere below that but i could be wrong.

h = 0.2;
k = 0.02;
alpha = 1;
lamda = alpha^2*k/h^2;
x = [0:h:2]';
t = [0:k:0.5]';
u = zeros(26,11);
u(:,1)=2;
u(:,11)=2;
for i=2:10
    u(1,i)=2*(x(i)-1).^3;
end
A = zeros(9,9);
B = zeros(9,9);
for i=1:9
A(i,i) = 2*(1+lamda);
B(i,i) = 2*(1-lamda);
if i<=8
A(i,i+1) = -lamda;
B(i,i+1) = lamda;
end
if i>=2
A(i,i-1) = -lamda;
B(i,i-1) = lamda;
end
end
for j=2:26
uu = A\(B*[u(j-1,2); u(j-1,3); u(j-1,4); u(j-1,5); u(j-1,6); u(j-1,7); u(j-1,8); u(j-1,9); u(j-1,10)]);
for kk=1:length(uu)
u(j,kk+1) = uu(kk);
end
end

surf(x,t,u)

Subject: Crank-Nicolson method

From: Torsten Hennig

Date: 23 Aug, 2010 12:10:47

Message: 16 of 17

The linear system of equations for the CN-method reads

u_{1}^(n+1) = u_{1}^{n}
-r*u_{i-1}^{n+1}+(1+2*r)*u_{i}^{n+1}-r*u_{i+1}^{n+1}=
r*u_{i-1}^{n}+(1-2*r)*u_{i}^{n}+r*u_{i+1}^{n}
(2<=i<=10)
u_{11}^{n+1} = u_{11}^{n}

with r = alpha*dt/(2*(dx)^2).

I can't find the corresponding coefficients
in the code you listed.

Best wishes
Torsten.

Subject: Crank-Nicolson method

From: Bill

Date: 23 Aug, 2010 21:55:55

Message: 17 of 17

On Aug 23, 2:52 pm, "Scott " <boofheadf...@hotmail.com> wrote:
> I took away the header function heat_crank and made a new m-file called f_x with the function in it and ran heat_crank and it appears to have worked. I have to show how the temperature (z-axis) appears over time and length so I want to make a mesh grid of the values of temperature over the time 0<t<0.5 and length 0<x<2. I can't get this to work. I'm not sure if all the temperature values have been calculated or if it's just calculated temperature for a certain x or t value...
>
> Any suggestions?

Hi Scott,

Did you end up getting this code to work correctly?

Regards,
Bill Hutton

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us